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Zusammenfassung 


In dieser Dissertation fiihre ich eine neue Observable der grossraumigen Struktur des 
Universums ein, das ortsabhangige Leistungsspektrum. Diese GroBe bietet ein Mass fur 
den “gequetschten” Limes der Dreipunktfunktion (Bispektrum), das heisst, eine Wellen- 
zahl ist wesentlich kleiner als die beiden iibrigen. Physikalisch beschreibt dieser Limes 
der Dreipunktfunktion die Modulation des Leistungsspektrums auf kleinen Skalen durch 
grossraumige Moden. Diese Modulation wird sowohl durch die Schwerkraft bewirkt als 
auch (moglicherweise) durch die kosmologische Inflation im friihen Universum. 

Fur die Messung teilen wir das Gesamtvolumen der Himmelsdurchmusterung, oder 
kosmologischen Simulation, in Teilvolumina ein. In jedem Teilvolumen messen wir die 
Uberdichte relativ zur mittleren Dichte der Materie (oder Anzahldichte der Galaxien) 
und das lokale Leistungsspektrum. Anschliessend messen wir die Korrelation zwischen 
Uberdichte und Leistungsspektrum. Ich zeige, dass diese Korrelation einem Integral fiber 
die Dreipunktfunktion entspricht. Wenn die Skala, an der das Leistungsspektrum aus- 
gewertet wird (inverse Wellenzahl, urn genau zu sein), vicl kleiner als die GroBe des Teil- 
volumens ist, dann ist das Integral fiber die Dreipunktfunktion vom gequetschten Limes 
dominiert. 

Lhn physikalisch zu verstehen, wie eine grossraumige Dichtefluktuation das lokale Leis- 
tungsspektrum beeinflusst, wenden wir das Bild vom “unabhangigen Universum” (“sep¬ 
arate universe”) an. Im Kontext der allgemeinen Relativitatstheorie kann eine lang- 
wcllige Dichtefluktuation exakt durch eine Friedmann-Robertson-Walker-(FRW-)Raumzeit 
beschrieben werden, deren Parameter sich von der “wahren” FRW-Raumzeit unterscheiden 
und eindeutig von der Dichtefluktuation bestimmt werden. Die Modulation des lokalen 
Leistungsspektrums kann dann durch die Strukturbildung innerhalb der modifizierten 
FRW-Raumzeit beschrieben werden. Insbesondere zeige ich, dass die Dreipunktfunktion 
im gequetschten Limes durch cliesen Ansatz einfacher und besser beschrieben wird als durch 
die herkommliche Herangehensweise mittels Storungstheorie. 

Diese neue Observable ist nicht nur einfach zu interpretieren (sie stellt die Antwort 
des lokalen Leistungsspektrums auf eine grofiskalige Dichtestorung dar), sie ermoglicht 
zudem die komplexe Berechnung der vollcn Dreipunktsfunktion zu umgehen, weil das Leis- 
tungsspektrum genauso wie die mittlere Dichte wesentlich lcichter als die Dreipunktsfunk¬ 
tion zu bestimmen sind. 

Anschliefiend wende ich die glciche Methodik auf die Daten der Himmelsdurchmus- 
tering SDSS-III Baryon Oscillation Spectroscopic Survey (BOSS) an, insbesondere den 
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Data Release 10 CMASS Galaxienkatalog. Wie ich zeige, stimmt das in den wirklichen 
Daten gemessene ortsabhangige Leistungsspektrum mit den sogenannten “mock” (also 
simulierten) Galaxienkatalogen iiberein, die auf dem PTHalo-Algorithmus basieren und 
die raumlichde Verteilung der wirklichen Galaxien im statistischen Sinne moglichst genau 
beschreiben wollcn. Genauer gesagt, liegen die Daten innerhalb der Streuung, die das 
ortsabhangige Leistungsspektrum zwischen den verschiedenen Realisicrungen von “mock” 
Katalogen aufweist. Diese Streuung betragt ca. 10% des Mittelewerts. In Kombina- 
tion mit dem (anisotropen) globalen Leistungsspektrum der Galaxien sowie dem Signal 
im schwachen Gravitationslinseneffekt, benutze ich diese 10%-Messung des ortsabhangigen 
Leistungsspektrums, urn den quadratischen Bias-Parameter der von BOSS gemessenen 
Galaxien zu bestinmicn, mit dem Ergebnis 62 = 0.41 ± 0.41 ( 68 % Vertrauensintervall). 

Schliefilich verallgemcinern wir die Analyse der Antwort des lokalen Leistungsspektrums 
auf eine Haufung von m grofiraumigen Wellenlangenmoden, wobei m < 3. In Analogie 
zum vorherigen Fall, kann die resultierende Modulation des Leistungsspektrums mit der 
m + 2-Punktskorrelationsfunktion im Limes gequetschter Konhgurationen (so dass immer 
zwei Wellenlangen wesentlich langer sind als die anderen), gemittled iiber die auftretenden 
Winkel, in Verbindung gebracht werden. Mit Hilfc von Simulationen “unabhangiger Uni- 
versen”, dass heiBt A-body-Simulationen in Anwesenhcit von Dichtestorungen unendlicher 
Lange, vergleichen wir unsere semianalytischen Modelle, die auf dem Bild der unabhangigen 
LIniversen basieren, mit den vollstandig nichtlinearen Simulationen bei bisher unerreichter 
Genauigkeit. Zudem testen wir die Annahme der gewohnlichen Storungstheorie, dass 
die nichtlineare N-Punktskorrelationsfunktion vollstandig durch das lineare Leistungsspek¬ 
trum bestimmt ist. Wir hnden bereits Abweichungen von 10% bei Wellenzahlen von 
k ~ 0.2 — 0.5 h Mpc -1 fur die Drei- bis Funf-Punktskorrelationsfunktion bei Rotver- 
schiebung z = 0. Dieses Ergebnis deutet darauf hin, dass die gewohnlichye Storungstheorie 
nicht ausreicht urn die Dynamik kollissionsloser Teilchen fur Wellenzahlen grdBer als diese 
korrekt vorherzusagen, selbst wenn alle hoheren Ordnungen in die Berechnung mit einbe- 
zogen werden. 



Abstract 


We present a new observable, position-dependent power spectrum, to measure the large- 
scale structure bispectrum in the so-called squeezed configuration, where one wavenumber, 
say k 3 , is much smaller than the other two, i.e. k 3 k\ ~ k 2 . The squeezed-limit bis- 
pectrum measures how the small-scale power spectrum is modulated by a long-wavelength 
scalar density fluctuation, and this modulation is due to gravitational evolution and also 
possibly due to the inflationary physics. 

We divide a survey volume into smaller subvolumes. We compute the local power 
spectrum and the mean overdensity in each smaller subvolume, and then measure the cor¬ 
relation between these two quantities. We show that this correlation measures the integral 
of the bispectrum, which is dominated by the squeezed configurations if the wavenumber 
of the local power spectrum is much larger than the corresponding wavenumber of the size 
of the subvolumes. This integrated bispectrum measures how small-scale power spectrum 
responds to a long-wavelength mode. 

To understand theoretically how the small-scale power spectrum is affected by a long- 
wavelength overdensity gravitationally, we use the “separate universe picture.” A long- 
wavelength overdensity compared to the scale of interest can be absorbed into the change 
of the background cosmology, and then the small-scale structure formation evolves in this 
modified cosmology. We show that this approach models nonlinearity in the bispectrum 
better than the traditional approach based on the perturbation theory. 

Not only this new observable is straightforward to interpret (the response of the small- 
scale power spectrum to a long-wavelength overdensity), but it also sidesteps the complexity 
of the full bispectrum estimation because both power spectrum and mean overdensity are 
significantly easier to estimate than the full bispectrum. 

We report on the first measurement of the bispectrum with the position-dependent 
correlation function from the SDSS-111 Baryon Oscillation Spectroscopic Survey (BOSS) 
Data Release 10 CMASS sample. We detect the amplitude of the bispectrum of the BOSS 
CMASS galaxies at 7.4<r, and constrain their nonlinear bias to be b 2 = 0.41 ± 0.41 (68% 
C.L.) combining our bispectrum measurement with the anisotropic clustering and the weak 
lensing signal. 

We finally generalize the study to the response of the small-scale power spectrum to m 
long-wavelength overdensities for m < 3. Similarly, this response can be connected to the 
angle-average (m+2)-point function in the squeezed configurations where two wavenumbers 
are much larger than the other ones. Using separate universe simulations, i.e. A-body sim- 
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ulations performed in the presence of an infinitely long-wavelength overdensity, we compare 
onr semi-analytical models based on the separate universe approach to the fully nonlinear 
simulations to unprecedented accuracy. We also test the standard perturbation theory 
hypothesis that the nonlinear n-point function is completely predicted by the linear power 
spectrum at the same time. We find discrepancies of 10% at k ~ 0.2 — 0.5 h Mpc -1 for 
three- to five-point functions at z — 0. This result suggests that the standard perturba¬ 
tion theory fails to describe the correct dynamics of collisionless particles beyond these 
wavenumbers, even if it is calculated to all orders in perturbations. 



Chapter 1 
Introduction 


1.1 Why study the position-dependent power spec¬ 
trum of the large-scale structure? 

The standard cosmological paradigm has been well developed and tested by the obser¬ 
vations of the cosmic microwave background (CMB) and the large-scale structure. The 
inhomogeneities seen in the universe originate from quantum fluctuations in the early uni¬ 
verse, and these quantum fluctuations were stretched to macroscopic scales larger than the 
horizon during the cosmic inflation [69j 11321 |5j [99], which is the early phase with expo¬ 
nential growth of the scale factor. After the cosmic inflation, the hot Big-Bang universe 
expanded and cooled down, and the macroscopic inhomogeneities entered into horizon and 
seeded all the structures we observe today. 

With the success of connecting the quantum fluctuations in the early universe to the 
structures we see today, the big questions yet remain: What is the physics behind inflation? 
Also, what is nature of dark energy, which causes the accelerated expansion in the late¬ 
time universe (see m for a review)? As the standard cosmological paradigm passes almost 
all the tests from the current observables, especially the two-point statistics of CMB and 
galaxy surveys, it is necessary to go to higher order statistics to obtain more information 
and critically test the current model. In particular, the mode coupling between a long- 
wavelength scalar density fluctuation and the small-scale structure formation receives much 
attention in the past few years. This coupling is due to the nonlinear gravitational evolution 
(see Ha for a review), and possibly the inflationary physics. Therefore, this provides 
a wonderful opportunity to test our understanding of gravity, as well as to probe the 
properties of inflation. 

Traditionally, the n-point function with n > 2 is used to characterize the mode cou¬ 
pling. Specifically, if one is interested in the coupling between one long-wavelength mode 
and two short-wavelength modes, we measure the three-point correlation function or its 
Fourier counterpart, the bispectrum, in the so-called “squeezed configurations,” in which 
one wavenumber, say k 3 , is much smaller than the other two, i.e. k$ -C k\ & k 2 - 

In the simplest model for the primordial non-Gaussianity (see [57] for a review on the 
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general primordial non-Gaussianities from various inflation models), the primordial scalar 
potential is given by 

$(r) = <j)( r) + / NL [0 2 (r) - (/ 2 (r))] , (1.1) 

where 0(r) is a Gaussian field and /nl is a constant characterizing the amplitude of the non- 
Gaussianity, which encodes the properties of inflation. Note that (0 2 (r)) assures (<f>(r)) = 
0. This simple model is known as the local-type primordial non-Gaussianity because <F(r) 
depends locally on /(r). The bispectrum of this local model is 


-£<&(ki, k 2 , k 3 ) = 2/ NL [P< s ,(k 1 )P< s ,(k 2 ) + 2 cyclic] , (1.2) 


where Pq>{k) oc k na ~ 4 is the power spectrum of the primordial scalar potential and n s ~ 0.96 
is its spectral index [88]. We can rewrite Bq> of this local model by fixing one wavenumber, 
say k \, as 


S$(ki,k 2 ,k 3 ) oc fc 2(ns 4) 


oc 


^2(n s 4) 






n s —4 


(1.3) 


where k 3 = —k 3 — k 2 because of the assumption of homogeneity. Bq> apparently peaks at 
k\ & k 2 fc 3 ~ 0, so the local-type primordial non-Gaussianity is the most prominent in 
the squeezed-limit bispectrum. 

Constraining the physics of inflation using the squeezed-limit bispectrum of CMB is a 
solved problem With the Planck satellite, the current constraint on the local-type 
primordial non-Gaussianity is /nl = 2.5 ±5.7 (68 % C.L.) using the temperature data alone 
and /nl = 0.8 ± 5.0 using the temperature and polarization data |128j . These are close to 
the best limits obtainable from CMB. To improve upon them, we must go beyond CMB to 
the large-scale structure, where observations are done in three-dimensional space (unlike 
CMB embedded on a two-dimensional sphere). Thus, in principle, the large-scale structure 
contains more information to improve the constraint on the physics of inflation. 

Measuring the three-point function from the large-scale structure (e.g. distribution of 
galaxies), however, is considerably more challenging compared to CMB. From the mea¬ 
surement side, the three-point function measurements are computationally expensive. In 
configuration space, the measurements rely on hireling particle triplets with the naive al¬ 
gorithm scaling as Ad) ar where IVp ar is the number of particles. Current galaxy redshift 
surveys contain roughly a million galaxies, and we need 50 times as many random sam¬ 
ples as the galaxies for characterizing the survey window function accurately. Similarly, in 
Fourier space, the bispectrum measurements require counting all possible triangle configu¬ 
rations formed by different Fourier modes, which is also computationally expensive. From 
the modeling side, galaxy surveys have more complicated survey selection function, which 
can bias the estimation (see e.g. |30j). Additionally, the nonlinear gravitational evolution 
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overdensity 


underdensity 


Figure 1.1: In the absence of the correlation between long and short wavelength modes, 
the blue fluctuation have the same statistical property in the overdense and underdense 
regions. On the other hand, in the presence of the positive correlation between long and 
short wavelength modes, the red fluctuation has larger (smaller) variance in the overdense 
(undersense) region. 


of matter density field and the complexity of galaxy formation make it challenging to ex¬ 
tract the primordial signal. The above difficulties explain why only few measurements of 
the three-point function of the large-scale structure have been reported in the literature 

Since our main interest is to measure the three-point function of the squeezed config¬ 
urations, there is a simpler way to sidestep all the above complexities of the three-point 
function estimation. As we stated, in the presence of the mode coupling between long- 
and short-wavelength modes, a long-wavelength density fluctuation modifies the small- 


sketches the short-wavelength modes with (red) and without (blue) correlation with the 
long-wavelength mode. As a consequence, for example, the n-point statistics and the halo 
mass function would depend on the local long-wavelength overdensity, or equivalently the 
position in space. Measurements of spatially-varying observables capture the effect of mode 
coupling, and can be used to test our understanding of gravity and the physics of infla¬ 
tion. A similar idea of measuring the shift of the peak position of the baryonic acoustic 
oscillation in different environments has been studied in m- 

In this dissertation, we focus on the position-dependent two-point statistics (see [321 
1115] for the mass function). Consider a galaxy redshift survey or simulation. Instead 
of measuring the power spectrum within the entire volume, we divide the volume into 
many subvolumes, within which we measure the power spectrum. These power spectra of 
subvolumes vary spatially, and the variation is correlated with the mean overdensities of 
the subvolumes with respect to the entire volume. This correlation measures an integral of 
the bispectrum (or the three-point function in configuration space), which represents the 


scale structure formation, and so the observables become position-dependent. Figure 1.1 
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response of the small-scale clustering of galaxies (as measured by the position-dependent 
power spectrum) to the long-wavelength density perturbation (as measured by the mean 
overdensity of the subvolumes). 

Not only is this new observable, position-dependent power spectrum , of the large-scale 
structure conceptually straightforward to interpret, but it is also simpler to measure than 
the full bispectrum, as the machineries for the two-point statistics estimation are well de¬ 
veloped (see [56] for power spectrum and [91] for two-point function) and the measurement 
of the overdensity is simple. In particular, the computational requirement is largely allevi¬ 
ated because we explore a subset of the three-point function corresponding to the squeezed 
configurations. More precisely, in Fourier space we only need to measure the power spec¬ 
trum, and in configuration space the algorithm of measuring the two-point function by 
finding particle pairs scales as A^ s (A^ par /iV s ) 2 = N^ ar /N s for the entire volume with N s 
being the number of subvolumes. In addition, for a fixed size of the subvolume, the mea¬ 
surement depends on only one wavenumber or one separation, so estimating the covariance 
matrix is easier than that of the full bispectrum from a realistic number of mock catalogs. 
The position-dependent power spectrum can thus be regarded as a useful compression of 
information of the squeezed-limit bispectrum. 

As this new observable uses basically the existing and routinely applied machineries to 
measure the two-point statistics, one can easily gain extra information of the three-point 
function, which is sensitive to the nonlinear bias of the observed tracers, from the current 
spectroscopic galaxy surveys. Especially, since the position-dependent power spectrum 
picks up the signal of the squeezed-limit bispectrum, it is sensitive to the primordial non- 
Gaussianity of the local type. 

This dissertation is organized as follows. In the rest of this chapter, we review the status 
of the observations of the large-scale structure, and the current theoretical understanding. 

In chapter [2j we introduce the main topic of this dissertation: position-dependent 
power spectrum and correlation function. We show how the correlation between the 
position-dependent two-point statistics and the long-wavelength overdensity is related to 
the three-point statistics. We also make theoretical template for this correlation using 
various bispectrum models. 

In chapter [3j we introduce the “separate universe approach,” in which a long-wavelength 
overdensity is absorbed into the background, and the small-scale structure formation 
evolves in the corresponding modified cosmology. This is the basis for modeling the 
response of the small-scale structure formation to the long-wavelength overdensity. We 
consider the fiducial cosmology to be flat ACDM, and show that the overdensity acts as 
the curvature in the separate universe. 

In chapter [4], we measure the position-dependent power spectrum from cosmological 
iV-body simulations. We compare various theoretical approaches to modeling the mea¬ 
surements from simulations, particularly the separate universe approach when the scales of 
the position-dependent power spectrum are much smaller than that of the long-wavelength 
overdensity. We also study the dependences of the position-dependent power spectrum on 
the cosmological parameters, as well as using the Fisher matrix to predict the expected 
constraints on biases and local-type primordial non-Gaussianity for current and future 
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galaxy surveys. 

In chapter [5j we report on the first measurement of the three-point function with the 
position-dependent correlation function from the SDSS-III Baryon Oscillation Spectro¬ 
scopic Survey Data Release 10 (BOSS DR10) CMASS sample. We detect the amplitude 
of the three-point function of the BOSS CMASS galaxies at 7.4cr. Combining the con¬ 
straints from position-dependent correlation function, global two-point function, and the 
weak lensing signal, we determine the quadratic (nonlinear) bias of BOSS CMASS galaxies. 

In chapter [6j we generalize the study to the response of the small-scale power spectrum 
in the presence of m long-wavelength modes for m < 3. This response can be linked to 
the angular-averaged squeezed limit of (m + 2)-point functions. We shall also introduce 
the separate universe simulations, in which Wbody simulations are performed in the pres¬ 
ence of a long-wavelength overdensity by modifying the cosmological parameters. The 
separate universe simulations allow unprecedented measurements for the squeezed-limit 
n-point function. Finally, we test the standard perturbation theory hypothesis that the 
nonlinear n-point function is completely predicted by the linear power spectrum at the 
same time. We find discrepancies of 10% at k ~ 0.2 — 0.5 h Mpc -1 for five- to three-point 
functions at z = 0. This suggests the breakdown of the standard perturbation theory, and 
quantifies the scales that the corrections of the effective fluid become important for the 
higher order statistics. 

In chapter [7] we summarize this dissertation, and present the outlook. 


1.2 Observations and measurements of the large-scale 
structure 

In the 1960s, before the invention of the automatic plate measuring machine and the 
densitometer, the galaxy catalogs such as Zwicky m and Lick [ 145] relied on visual 
inspection of poorly calibrated photographic plates. These surveys consisted of different 
neighboring photographic plates, so the uniformity of the calibration, which might cause 
large-scale gradients in the observed area, was a serious issue. Because of the lack of the 
redshifts (depths) of galaxies, only the angular clustering studies were possible. In addition, 
the sizes of the surveys were much smaller than the ones today, thus only the clustering 
on small scales, where the nonlinear effect is strong, can be studied. Nevertheless, in the 
1970s Peebles and his collaborators did the first systematic study on galaxy clustering using 
the catalogs at that time. The series of studies, starting with |125j . considered galaxies 
as the tracers of the large-scale structure for the first time, which was a ground-breaking 
idea. These measurements confirmed the power-law behavior of the angular two-point 
function, and the interpretation was done in the framework of Einstein-de Sitter universe, 
i.e. matter-dominated flat Friedmann-Lemaitre-Robertson-Walker universe. 

In the 1980s, the invention of the automatic scanning machines as well as CCDs revolu¬ 
tionized the large-scale structure surveys, and resulted in a generation of wide-held surveys 
with better calibration and a three-dimensional view of the universe. Photographic plates 
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Figure 1.2: (Credit: Michael Blanton and SDSS collaboration) SDSS galaxy maps projected 
on the RA-redshift plane. (Left) The SDSS main galaxy sample out to z ~ 0.15. (Right) 
Yellow, red, and white dots are the SDSS main galaxy sample, luminous red galaxies, and 
CMASS sample, respectively, out to z ~ 0.7. 


became obsolete for the large-scale structure studies, and nowadays photometric surveys 
use large CCD cameras with millions of pixels. The galaxy redshift surveys, which gener¬ 
ally require target selections with photometric detection and then spectroscopic follow-up, 
thus open a novel avenue to study the universe. It was shown that the redshift-space 
two-point correlation function in the CfA survey ED agreed well with the previous studies 
on angular clustering, if the redshift direction is integrated over [46]. In the 1990s, the 
number of galaxies in surveys was ~ 10 3 — 10 4 , but with these data it was already shown 
that the large-scale power spectrum was inconsistent with the CDM model [52l 11351 168] . 
in agreement with the study done in the angular clustering HQU- 


Another quantum leap of the sizes of the galaxy surveys happened in the 2000s, when 
the technology of the massive multi-fiber or multi-slit spectroscopy became feasible. Sur¬ 
veys such as Two-degree-Ficld Galaxy Redshift Survey (2dFGRS) [34] and Sloan Digital 
Sky Survey (SDSS) [ 172] targeted at obtaining spectra of ~ 10 5 — 10 6 galaxies. The left 
panel of figure L2 shows the SDSS main galaxy sample out to z ~ 0.15, which corresponds 
to roughly 440 h^ 1 Mpc. It is clear even visually that the distribution of galaxies follows 
filamentary structures, with voids in between the filaments. These data contain precious 
information of the properties of the universe. For example, in 2005, the baryonic acoustic 
oscillations (BAO) in the two-point statistics were detected for the first time by 2dFGRS 
in the power spectrum and SDSS in the two-point correlation function, as shown 
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Figure 1.3: (Credit: 2dFGRS and SDSS collaborations) The first detection of BAOs in 
2dFGRS (left) and SDSS (right). The analysis of 2dFGRS was done in Fourier space, and 
the BAOs are the wiggles in the power spectrum at 0.05 h Mpc” 1 < k < 0.2 h Mpc” 1 ; 
the analysis of SDSS was done in configuration space, and the BAOs are the bump in the 
two-point correlation function at r r\j 100 h 1 Mpc. 


in figure L3 The detection is phenomenal given the fact that the BAOs are of order a few 
percent features on the smooth functions. Galaxy surveys in this era started probing the 
weakly nonlinear regime, where the theoretical understanding is better, so we can extract 
cosmological information. 

With the success of the first BAO detection, more galaxy redshift surveys, such as 
WiggleZ [22] [23] and SDSS-III Baryon Oscillation Spectroscopic Survey (BOSS) 0 [6], 


followed and extended to higher redshifts. The right panel of figure 1.2 shows the galaxy 
map of SDSS out to z ~ 0.7, which corresponds to roughly 1800 h~ l Mpc. The BAO 
feature can be used as a standard ruler to measure the angular diameter distance and 
Hubble expansion rate. This is particularly useful for studying the time-dependence of 
dark energy, which began to dominate the universe at z ~ 0.4 where the galaxy clustering 
is measured. The BAOs in the galaxy clustering thus becomes a powerful probe of dark 
energy. 

Thus far, most of the studies have focused on the two-point statistics. However, there 
is much more information in the higher-order statistics. Especially, ongoing galaxy surveys 
such as Dark Energy Survey [162J and the extended BOSS, as well as upcoming galaxy 
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surveys such as Hobby Eberly Telescope Dark Energy experiment (HETDEX) [Tlj and 
Subaru Prime Focus Spectrograph [ 156 ] will measure the galaxies at even higher redshift. 
For instance, HETDEX will use Lyman-alpha emitters as tracers to probe the matter 
distribution at 1.9 < z < 3.5. At such high redshift, the gravitational evolution is relatively 
weak and can still be predicted analytically, hence this is an ideal regime to critically test 
our understanding of gravity, as well as the physics of inflation via the primordial non- 
Gaussianity. 

Currently, most of the constraint on the local-type primordial non-Gaussianity from 
the large-scale structure is through the scale-dependent bias m HOG H51j . That is, dark 
matter halos (or galaxies) are biased tracers of the underlying matter distribution forming 
at density peaks, so the formation of halos would be modulated by the additional correlation 
between the long and short wavelength modes due to the primordial non-Gaussianity. As 
a result, the halo bias contains a k~ 2 scale-dependent correction, and this correction is 
prominent on large scale. Measurements of the large-scale galaxy power spectrum, which is 
proportional to bias squared, can thus be used to constrain the primordial non-Gaussianity 
of the local type. 

This distinct feature appears in the galaxy power spectrum at very large scales, hence it 
is crucial to have a huge survey volume to beat down the cosmic variance. Moreover, if the 
galaxies are highly biased, then the signal-to-noise ratio would also increase. Thus, many 
studies have used quasars at 0.5 < z < 3.5 from SDSS to constrain the primordial non- 
Gaussianity mmm- Similar methods, such as combining the abundances and clustering 
of the galaxy clusters cna as well as the correlation between CMB lensing and large-scale 
structure [EIIISO], have also been proposed to study the primordial non-Gaussianity from 
large-scale structure. It is predicted in [92] that for Large Synoptic Survey Telescope [100] 
the constraint on /nl, parametrization of the local-type primordial non-Gaussianity, using 
the scale-dependent bias can reach ct(/nl) ~ 5 (95% C.L.). 

The error bar on /nl from the scale-dependent bias is limited by the number of Fourier 
modes on large scales. On the other hand, for the bispectrum analysis, we are looking 
for triangles formed by different Fourier modes, so the bispectrum contains more infor¬ 
mation and will have a tighter constraint on /n l- The difficulty for using the large-scale 
structure bispectrum to constrain / NL is that gravity produces non-zero squeezed-limit, bis¬ 
pectrum even without primordial non-Gaussianity, and the signal from gravity dominates 
for the current limit on /nl- This is why recent measurements of the large-scale structure 
bispectrum have focused on constraining the growth and galaxy biases mm- 

The difficulty in modeling nonlinear effects can be alleviated if the observations are done 
in the high-redshift universe, where the gravitational evolution on quasi-linear scales can 
still be described by the perturbation theory approach. While theoretically we are reaching 
the stage for studying the higher-order statistics, e.g. the three-point correlation function 
or the bispectrum, if the data are obtained at high redshift, in practice the measurements 
and analyses are still computational challenging. It is thus extremely useful to fold a way 
to compress the information, such that studying the three-point correlation function of the 
galaxy clustering is feasible. 

Another complexity of the bispectrum measurement from galaxy surveys is the window 
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function effect. Namely, galaxy surveys almost always have non-ideal survey geometry, 
e.g. masking around the close and bright objects or the irregular boundaries, as well as 
the spatial changes in the extinction, transparency, and seeing. These effects would bias 
the measurement, and extracting the true bispectrum signal becomes difficult. While the 
observational systematics would enters into the estimation of both power spectrum and 
bispectrum, the technique of deconvolving the window function effects, e.g. |134l 1133] . has 
been relatively well developed for the two-point statistics. 

The subject of this dissertation is to find a method to more easily extract the bispec¬ 
trum in the squeezed configurations, where one wavenumber, say Aq, is much smaller than 
the other two, i.e. Aq -C fci ~ Aq. The squeezed-limit bispectrum measures the correlation 
between one long-wavelength mode (Aq) and two short-wavelength modes (Aq, Aq), which is 
particularly sensitive to the local-type primordial non-Gaussianity. Specifically, we divide 
a survey into subvolumes, and measure the correlation between the position-dependent 
two-point statistics and the long-wavelength overdensity. This correlation measures an in¬ 
tegral on the bispectrum, and is dominated by the squeezed-limit signal if the wavenumber 
of the position-dependent two-point statistics is much larger than the wavenumber cor¬ 
responding to the size of the subvolumes. Therefore, without employing the three-point 
function estimator, we can extract the squeezed-limit bispectrum by the position-dependent 
two-point statistics technique. Furthermore, nonlinearity of the correlation between the 
position-dependent two-point statistics and a long-wavelength mode can be well modeled 
by the separate universe approach, in which the long-wavelength overdensity is absorbed 
into the background cosmology; the window function effect can also be well taken care of 
because this technique measures essentially the two-point function and the mean overden¬ 
sity, for which the procedures for removing the window function effects are relatively well 
developed. With the above advantages, the position-dependent two-point statistics is thus 
a novel and promising method to study the squeezed-limit bispectrum of the large-scale 
structure. 

Galaxy redshift surveys has entered a completely new era at which the sizes (e.g. survey 
volume and number of observed galaxies) are huge, and redshifts are high. As the signal-to- 
noise ratio of the higher-order statistics, especially for the primordial non-Gaussianity, will 
be much higher in the upcoming galaxy surveys than the previous ones, we should do our 
best to extract the precious signal for improving our understanding of the universe. The 
new observable, position-dependent power spectrum, proposed in the dissertation would 
help us achieve this goal. 

1.3 Theoretical understanding of the large-scale struc¬ 
ture 

1.3.1 Simulations 

How do we understand the gravitational evolution of the large-scale structure? Because 
of the process is nonlinear, the gold standard is the cosmological iV-body simulations of 
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Figure 1.4: (Credit: Volker Springel) A slice of the Millennium simulation shows the cosmic 
web of the large-scale structure. 


collisionless particles. 

Using IV-body simulations to solve for gravitational dynamics has a long history (see 
e.g. [4B] for a review). The first computer calculation was done back in 1963 by [lj 
with N = 16. Later, N has roughly doubled every two years following Moore’s law, and 
nowadays state-of-the-art simulations for collisionlcss particles have N = 10 9 — 5 x 10 10 , 
e.g. |154l 11611 1551 jgBlj . The IV-body codes, such as GADGET-2 [T53j . solve dynamics 
of dark matter particles, and dark matter particles are grouped into dark matter halos 
by algorithms such as friends-of-friends or spherical overdensity. We can thus study the 
properties of these halos, such as the clustering and the mass function. Figure |L4| shows 
a slice of the Millennium simulation [154] . The cosmic web of the large-scale structure is 
obvious, and we also find the similarity between the simulation and the observation, i.e. 


figure 1.2 


Direct observations, however, can only be done for luminous objects, e.g. galaxies, 
so we must relate galaxies to the halos. As baryonic physics is more complicated than 
gravity-only dynamics, simulations of galaxy formation can only be done with the usage 
of sub-grid physics models. That is, the empirical relations of the feedback from baryonic 
physic such as supernovae and active galactic nuclei are used in the galaxy formation sim¬ 
ulations. State-of-the-art hydrodynamic simulations for the galaxy formation are Illustris 
[ 169] and EAGLE [136] . Since these simulations are extremely computationally intensive, 
the simulation box size cannot be too large (e.g. ~ 100 h Mpc for both Illustris and 
EAGLE). On the other hand, galaxy redshift surveys at present day have sizes of order 
1-10 h~ 3 Gpc 3 . 


Alternatively, we can link the simulated dark matter halos to galaxies using techniques 
such as semi-analytic models | 81 j . which use the merging histories of dark matter halos, or 
halo occupation distribution [IB] ED], which uses the statistical relation between halos and 
galaxies. As both methods contain free parameters in the models, these parameters can be 
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tuned such that multiple properties (e.g. clustering and environmental dependence) of the 
“simulated” galaxies match the observed ones (see [68] for the recent comparison between 
semi-analytic models and the observations). These thus provide a practically feasible way 
to populate galaxies in dark matter halos in the simulations. 

The remaining task, especially for the clustering analysis of galaxy redshift surveys, 
is to generate a large suite of simulations with halos, so they can be used to estimate 
the covariance matrices of the correlation function or the power spectrum, which are the 
necessary ingredient for the statistical interpretation of the cosmological information. In 
particular, the present-day galaxy redshift surveys contain a huge volume, and if high 
enough mass resolution for halos (which depend on the properties of the observed galaxies) 
is required, iV-body simulation are not practical. Thus, most analyses for the two-point 
statistics of galaxy surveys use algorithms such that the scheme of solving dynamics is 
simplified to generate “mock” halos. For example, COLA [159j solves the long-wavelength 
modes analytically and short-wavelength modes by iV-body simulations, and PTHalos 
[13511TM [105] is based on the second-order Lagrangian perturbation theory. 

As the full TV-body simulations are regarded as the standard, we can input identical 
initial conditions to various codes for generating mock catalogs (halos) and compare the 
performances. Currently most of the observations now have been focused on the two-point 
statistics and the mass function, so these algorithms are designed to recover these two 
quantities. On the other hand, for the three-point function, which is more sensitive to 
nonlinear effects, more careful and systematic studies are necessary. One recent compari¬ 
son between major methodologies for generating mock catalogs shows that the differences 
between W-body simulations and mock generating codes are much larger for the three- 
point function than for the two-point function [31]. This suggests that at this moment the 
full W-body simulations are still required to understand the three-point function, or the 
bispectrum in Fourier space. 

The full bispectrum contains triangles formed by different Fourier modes. While the 
configurations of the bispectrum can be simulated easily if three wavenumbers are similar, 
the squeezed triangles are more diffi cult, because a large volume is needed to simulate the 
coupling between long-wavelength modes and small-scale structure formation. In partic¬ 
ular, if high enough mass resolution is required, the simulations become computationally 
demanding. In this dissertation, we provide a solution to this problem. Specifically, we 
absorb the long-wavelength overdensity into the modified background cosmology (which 
is the subject in chapter [3]) and perform the A r -body simulations in the separate universe 
(which is the subject in section 6.2). This setting simulates how the small-scale structure 
is affected by a long-wavelength mode. As the box size of the separate universe simula¬ 
tions can be small (~ 500 h~ 1 Mpc), increasing the mass resolution becomes feasible. This 
technique is therefore useful for understanding the nonlinear coupling between long and 
short wavelength modes. 

Another computational challenge to the bispectrum analysis is the number of Fourier 
bins. Specifically, the bispectrum contains all kinds of triangles, so the number of bins is 
much larger than that of the power spectrum. If the mock catalogs are used to estimate the 
covariance matrix, many more realizations are required to characterize the bispectrum than 
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the power spectrum. The lack of realizations of mock catalogs would result in errors in the 
covariance matrix estimation, and the parameter estimation would be affected accordingly 
[51] . Therefore, even if there is an algorithm to generate mock catalogs with accurate two- 
and three-point statistics, we still need a huge amount of them for data analysis, which can 
be computational challenging. The advantage of the position-dependent power spectrum 
is that it depends only on one wavenumber, so the number of bins is similar to that of 
the power spectrum. This means that we only need a reasonable number of realizations 
(~ 1000 ) for analyzing the power spectrum and position-dependent power spectrum jointly. 

1.3.2 Theory 

While iV-body simulations are the gold standard for understanding nonlinearity of the 
large-scale structure, it is impractical to run various simulations with different cosmological 
parameters or models. It is therefore equally important to develop analytical models so 
they are easier to compute. We can then use them in the cosmological inferences, e.g. the 
Markov chain Monte Carlo methods. 

The most commonly used technique is the perturbation theory approach, in which 
the fluctuations are assumed to be small so they can be solved recursively (see [T7] for a 
review). For example, in the standard perturbation theory (SPT), the density fluctuations 
and peculiar velocity fields are assumed to be small. Using this assumption, we can expand 
the continuity, Euler, and Poisson equations at different orders, and solve the coupled 
differential equations order by order (see appendix [A] and E3 for a brief overview on 
SPT). The SPT power spectrum at the first order contains the product of two first order 
fluctuations (Pn); the next-to-lcading order SPT power spectrum contains the products 
of Erst and third order fluctuations (P 13 ) as well as two second order fluctuations (P 22 )- A 
similar approach can be done in Lagrangian space, in which the displacement held mapping 
the initial (Lagrangian) position to the final (Eulerian) position of fluid element is solved 
perturbatively [ 1051, 107 ]. 

Generally, the perturbation theory approach works well on large scales and at high 
redshift, where nonlinearity is small. On small scale and at low redshift, the nonlinear effect 
becomes more prominent, including higher order corrections is thus necessary. However, 
the nonlinear effect can be so large that even including more corrections does not help. The 
renormalized perturbation theory (RPT) was introduced to alleviate the problem [TT[ ITS], 
Specifically, RPT categorizes the corrections into two kinds: the mode-coupling effects 
and the renormalization of the propagator (of the gravitational dynamics). Thus, in RPT 
the corrections for nonlinearity become better defined, and so the agreement with the 
nonlinear power spectrum extends to smaller scales compared to SPT. Another approach 
is the effective held theory (EFT) [25]: on large-scale the matter fluid is characterized by 
parameters such as sound speed and viscosity, and these parameters are determined by 
the small-scale physics that is described by the Boltzmann equation. In practice, these 
parameters are measured from ./V-body simulations with a chosen smoothing radius. As 
for RPT, EFT also gives better agreement with the nonlinear power spectrum on smaller 
scales compared to SPT. 
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A different kind of approach to compute the clustering properties of the large-scale 
structure is to use some phenomenological models. The well known phenomenological 
model is the halo model (see [35] for a review), where all matter is assumed to be contained 
inside halos, which are characterized by the density profile (e.g. NFW profile |118j ) and 
the mass function (e.g. Sheth-Tormen mass function [ 149] ). The matter 77 ,-point functions 
is then the sum of one-halo term (all n positions are in one halo), two-halo term ((77 — 1) 
positions are in one halo and the other one is in a different halo), to 77 ,-halo term (77 positions 
are in n different halos). The small-scale nonlinear matter power spectrum should thereby 
be described by the halo properties, i.e. the one-halo term. Some recent work attempted 
to extend the halo model to better describe the nonlinear matter power spectrum. For 
example, in [ 116] , the Zeldovich approximation [ 173] is added in the two-halo term, and 
a polynomial function (A 0 + A 2 k 2 + A 4 fc 4 + • • •) is added to model the baryonic effects; 
in [145] . the two-halo term with the Zeldovich approximation is connected to the SPT 
one-loop power spectrum (P u + P 13 + P 22 ), but a more more complicated function is used 
for the one-halo term. 

One can also construct fitting functions based on results of IV-body simulations. The 
most famous fitting functions of the nonlinear matter power spectrum are the halofit pre¬ 
scription [152] and the Coyote emulator [70]. More specifically, the Coyote emulator was 
constructed with a suite of IV-body simulations with a chosen range of cosmological pa¬ 
rameters (e.g. fib, h2 m ). Then, the power spectrum is computed based on the interpolation 
of the input cosmological parameters. Thus, the apparent limitation of these simulation- 
calibrated fitting formulae is that they are only reliable within limited range of cosmological 
parameters and restricted cosmological models. 

Similar to that of simulations, most of the theory work has been focused on precise 
description of nonlinearity of the matter power spectrum, while relatively few work has been 
devoted to investigate nonlinearity of the bispectrum. Therefore, the matter bispectrum 
is normally computed at the SPT tree-level, i.e. the product of two first order fluctuations 
and a second order fluctuation. Some studies pmiMifia] considered nonlinearity of the 
bispectrum by replacing the SPT kernel with fitting formulae containing some parameters, 
which are then obtained by fitting to IV-body simulations. These models, however, lack 
the theoretical foundation, and it is unclear how the fitting parameters would depend on 
the cosmological parameters. 

In this dissertation, we provide a semi-analytical model for describing the bispectrum 
in the squeezed configurations. Specifically, we show that the real-space angle-average 
squeezed-limit bispectrum is the response of the power spectrum to an isotropically in¬ 
finitely long-wavelength overdensity. Due to the presence of this overdensity, the back¬ 
ground cosmology is modified, and the small-scale power spectrum evolves as if matter is 
in the separate universe. In chapter [4| we show it is straightforward to combine the separate 
universe approach and the power spectrum computed from perturbation theory approach, 
phenomenological models, or simulation-calibrated fitting formulae. More importantly, the 
results of the separate universe approach agree better with the IV-body simulation mea¬ 
surements in the squeezed limit than that of the real-space bispectrum fitting formula. In 
chapter [6j we generalize the separate universe approach to the response of the small-scale 
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power spectrum to m infinitely long-wavelength overdensities for m < 3. As expected, this 
response is related to the squeezed-limit (m + 2)-point function with a specific configuration 
shown in figure 6.1 The separate universe approach is thus extremely useful for modeling 


the squeezed-limit n-point functions, and its analytical form can be added to the fitting 
formulae. 

As future galaxy surveys contain data with unprecedented amount and quality which 
can be used to test our understanding of gravity and the physics of inflation, accurate 
theoretical model, especially for the bispectrum, is required to achieve the goal. While 
the full bispectrum contains various triangles formed by different Fourier modes, in this 
dissertation we present the theoretical model specifically for the squeezed triangles, so more 
work needs to be done for the other configurations. 




Chapter 2 

Position-dependent two-point 
statistics 


2.1 In Fourier space 


2.1.1 Position-dependent power spectrum 

Consider a density fluctuation field, <5(r), in a survey (or simulation) of volume V r . The 
mean overdensity of this volume vanishes by construction, i.e. 

5 — f d 3 r <5(r) = 0 . (2.1) 

V r JVr 

The global power spectrum of this volume can be estimated as 

Pfk) = T|i(k)| 2 , (2.2) 


where <5(k) is the Fourier transform of <5(r). 

We now identify a subvolume Vl centered at r^. The mean overdensity of this subvol¬ 
ume is 


<5(i~l) = tt [ d 3 r 6(r) 
Vl Jv l 


1 

Vl 


d 3 r <5(r)W(r — r L ) , 


(2.3) 


where W{ r) is the window function. For simplicity and a straightforward application to 
the IV-body simulation box, throughout this dissertation we use a cubic window function 
given by 


3 

W( r) = W L ( r) = 9(n), e(n) 

1=1 


l, \n\ < L/ 2, 
0 , otherwise . 


(2.4) 


where L is the side length of Vl- The results are not sensitive to the exact choice of 
the window function, provided that the scale of interest is much smaller than L. While 
<5 = 0, <5(r^) is non-zero in general. In other words, if 5(r^) is positive (negative), then this 
subvolume is overdense (underdense) with respect to the mean density in V r . 
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Similar to the definition of the global power spectrum in V r , we define the position- 
dependent power spectrum in Vl as 

-P(k, r L ) = -j-|<5(k, r L )| 2 , (2.5) 

Vl 

where 

<5(k, r L ) = f d 3 r <5(r)e _ * r ' k (2.6) 

Jv L 

is the local Fourier transform of the density fluctuation field. The integral ranges over 
the subvolume centered at r^. With this quantity, the mean density perturbation in the 
subvolume centered at r L is given by 

<5(i-l) = ^-6{k = 0,r L ) . (2.7) 

One can use the window function Wl to extend the integration boundaries to infinity as 

5( k, r L ) = f d 3 r 5{v)W L {v - r L )e~ 4r k = J ^ 5( k - q)W L {q)e*™ , (2.8) 

where Wl{ q) = L 3 f([? =1 sinc(gjL/2) is the Fourier transform of the window function and 
sinc(a;) = sin(x)/a:. Therefore, the posit ion-dependent power spectrum of the subvolume 
Vl centered at r L is 

P(k, r l) = y L / ^ ^ ~ qO'H-k - q2)tr i (q 1 )W' i (q 2 )e-‘'-(-‘ + ®> . (2.9) 


2.1.2 Integrated bispectrum 

The correlation between P(k,ri) and 5(r^) is given by 


(P(k,r L )5(r L )) 


1 r d 3 q 1 

VlJ (2vr) 3 


WfSwf < ' S(k_qi) ' 5( “ k " q2 )' 5 (“ q3 )> 

X W L ( qi )W L (q 2 )W L (q 3 )e- ir ^+v+^ , (2.10) 


where ( ) denotes the ensemble average over many universes. In the case of a simulation 
or an actual survey, the average is taken instead over all the subvolumes in the simulation 
or the survey volume. Through the definition of the bispectrum, (5(qi)<f(q2)^(q3)) = 
5(qi ! q2,q3)(2vr) 3 5 D (q 1 + q 2 + qs) where S D is the Dirac delta function, eq. ( |2.10[ ) can be 
rewritten as 


(P(k,r L )5(r L )> =—= 


d 3 qi f d?q 3 


VlJ ( 2 tt ) 3 ./ 


= iB l ( k) . 


P(k - qi, -k + qi + q 3 , -q 3 ) 

xkFL(qi)fk L (- qi -q3)kFL(q3) 


( 2 . 11 ) 
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As anticipated, the correlation of the position-dependent power spectrum and the local 
mean density perturbation is given by an integral of the bispectrum, and we will therefore 
refer to this quantity as the integrated bispectrum , iBi( k). 

As expected from homogeneity, the integrated bispectrum is independent of the location 
(r/,) of the subvolumes. Moreover, for an isotropic window function and bispectrum, the 
result is also independent of the direction of k. The cubic window function eq. (2.4) is 
of course not entirely spherically symmetric [j] and there is a residual dependence on k in 
eq. (2.11). In the following, we will focus on the angle average of eq. (2.11), 


iB L (k ) = 





^ P(k,r L ) ) 5(t l ) 


B( k - qi, -k + qi + q 3 , -q 3 ) 
xW L ( qi )W L (- qi -q 3 )W L (q 3 ) 


( 2 . 12 ) 


The integrated bispectrum contains integrals of three sine functions, sine (a;), which are 
damped oscillating functions and peak at \x\ < tt. Most of the contribution to the in¬ 
tegrated bispectrum thus comes from values of q\ and g 3 at approximately 1/L. If the 
wavenumber k we are interested in is much larger than 1/L (e.g., L = 300 hr 1 Mpc and 
k > 0.3 h Mpc -1 ), then the dominant contribution to the integrated bispectrum comes 
from the bispectrum in squeezed configurations, i.e., P(k — qi,—k + qi + q 3 , — q 3 ) —> 
B( k, —k, —q 3 ) with qi -C k and g 3 <C k. 


2.1.3 Linear response function 

Consider the following general separable bispectrum, 

B( kx, k 2 , k 3 ) = /(ki, k 2 )P(k 1 )P(k 2 ) + 2 cyclic, (2.13) 

where /(ki, k 2 ) = f(k\, fc 2 , k\ ■ k 2 ) is a dimensionless symmetric function of two k vectors 
and the angle between them. If / is non-singular as one of the k vectors goes to zero, we 
can write, to lowest order in q 1 /k and q^/k, 

B{ k - qi, -k + qi + q 3 , -q 3 ) = /(k - qi, -q 3 )P(|k - qi|)P(g 3 ) 

+ / (~ k + qi + q3, — qs)-P(| — k + q! + q3|)P(o , 3) 

+ /(k - qi, -k + qx + q 3 )P(|k - qi|)P(| - k + q 1 + q 3 |) 

= 2/(k, 0 )P(k)P(q 3 ) + /(k, —k)[P(fc)] 2 + O (/«) . 

(2.14) 

choose the cubic subvolumes merely for simplicity. In general one can use any shapes. For example, 
one may prefer to divide the subvolumes into spheres, which naturally lead to an isotropic integrated 
bispectrum iB^k). 
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2. Position-dependent two-point statistics 


For matter, momentum conservation requires that /(k, —k) = 0 [126], as can explicitly be 
verified for the P 2 kernel of perturbation theory. We then obtain 


d 2 k 


^ B( k - qi, -k + qi + q 3 , -q 3 ) = f(k)P(k)P(q 3 ) + O 


(2.15) 


where f[k) = 2/(0, k). Note that the terms linear in g 1)3 cancel after angular average. 
Since the window function in real space satisfies kF|(r) = Wl{ r), we have 


/ (^3 ~ q.s) = VF L (q 3 


Performing the q 3 integral in eq. (2.12) then yields 

Wl(q 3 )P(q 3 )f(k)P(k) = a 2 L f(k)P(k) , 


. / 7 \ kL —^oo 1 

iB L {k) = —= 


rf3( ?3 Tir2/ 


VI J (2vr)3 

where cp 2 l is the variance of the density field on the subvolume scale, 


cr? = 


„,2 


VIS (2jt)3 


Wl(q 3 )P(q 3 ) . 


(2.16) 


(2.17) 


(2.18) 


Eq. (2.17) shows that the integrated bispectrum measures how the small-scale power spec¬ 
trum, P(k), responds to a large-scale density fluctuation with variance cr/, with a response 
function given by f(k). 

An intuitive way to arrive at the same expression is to write the response of the small- 
scale power spectrum to a large-scale density fluctuation as 


P(k,r L ) = P(k)| 5=0 + 


dP(k) 


d5 


<*(ri) + • • • , 


(2.19) 


< 5=0 


where we have neglected gradients and higher derivatives of 5(r l). We then obtain, to 
leading order, 


iB L (k) = a 2 L 


2 d In P(k) 


dd 


P(k). 


( 2 . 20 ) 


< 5=0 


Comparing this result with eq. (2.17), we find that f{k) indeed corresponds to the lin¬ 
ear response of the small-scale power to the large-scale density fluctuation, din P(k)/S. 
Inspired by eq. (2.20), we define another quantity, the normalized integrated bispectrum, 

iB L {k) 


o\P(k) 


( 2 . 21 ) 


This quantity is equal to f(k) and the linear response function in the limit of kL —> oo. 
For the standard perturbation theory kernel 


/(ki, k 2 ) — P2(k l5 k 2 ) — - + - 


1 ki • k 2 / k 




W + + = 


7 2 k\k 2 \k 2 kj 7 V kik 2 


2 /ki • k 2 


2 


( 2 . 22 ) 
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in the squeezed limit the integrated bispectrum becomes (see appendix |A.2| for the detailed 
derivation) 


iB L {k) 


kL—toe 


47 

21 


1 din Pi(k) 
3 d In k 


Pi{k)a 


2 

L i 


(2.23) 


and the response function is 


/(*) 


47 _ 1 din Pi(k) 
21 3 dIn k 


We shall discuss more details in section 14.2.11 


(2.24) 


2.1.4 Integrated bispectrum of various bispectrum models 


To evaluate the integrated bispectrum, we insert the bispectrum models into eq. (2.12) and 


perform the eight-dimensional integral. Because of the high dimensionality of the integral, 
we use the Monte Carlo integration routine in GNU Scientific Library to numerically eval¬ 
uate iB L {k). Let us consider the simplest model of galaxy bispectrum with local-type 
primordial non-Gaussianity 


B g { ki, k 2 , k 3 ) — 6^5 S p T (k 1 , k 2 , k 3 ) + b\b 2 B b2 { k 1; 


k2, k 3 ) + &?/ N L-B/ NL (ki, k 2 , k 3 ) , (2.25) 


where b\ is the linear bias, b 2 is the quadratic nonlinear bias, and /nl is the parametrization 
for the local-type primordial non-Gaussianity. Note that the scale-dependent bias due to 
the local type non-Gaussianity (45^ 11061 1151 ] is neglected in eq. (2.25), and the latest 


bispectrum model with primordial non-Gaussianity can be found in [TO] 1158] . 


The first two terms of eq. (2.25) are due to the nonlinear gravitational evolution. 


Specifically, the standard perturbation theory (SPT) with local bias model predicts 


,see 


appendix A.2 for detailed derivation) 


5 SPT (k 1 , k 2 , k 3 ) = 2F 2 (ki, k 2 )Pi(k 1 , a)P(k 2 , a) + 2 cyclic 
B b2 ( ki, k 2 , k 3 ) = Pi(h, a)Pi(k 2 , a) + 2 cyclic , 


(2.26) 


where Pi(k) is the linear bispectrum. For the bispectrum of local-type primordial non- 
Gaussianity, we consider the local ansatz for the primordial scalar potential as [89] 

<L(r) = /(r) + /nl[</> 2 ( r) - (</> 2 (r))] , (2.27) 

where </>(r) is a Gaussian held, and /nl is a constant characterizing the amplitude of the 
primordial non-Gaussianity. As the density fluctuations are linked to the scalar potential 
through the Poisson equation 

i(k,o) = M(fc,o)*(k), M(k,a) = , (2.28) 

with D(a) and T[k ) being the linear growth factor and the transfer function respectively, 
in the leading order the primordial non-Gaussianity appears in matter bispectrum as 


f NL B fNL ( ki,k 2 ,k 3 ) = 2f NL M(k 1 ,a)M(k 2 ,a)M(k 3 ,a)[P^(k 1 )P^(k 2 ) + 2 cyclic] , (2.29) 
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2. Position-dependent two-point statistics 



Figure 2.1: Normalized integrated bispectrum in different sizes of subvolumes. The colored 
lines show non-Gaussian components at different redshifts assuming / NIj = 50, while the 
solid and dashed lines show linear and nonlinear bias components assuming bi = Z 2 = 1- 
The cut-off at low-fc corresponds to the fundamental frequencies of the subvolumes, 27 t/L. 


where P<f>(k) is the power spectrum of the scalar potential. 

Figure [27T| shows the normalized integrated bispectrum, for which we shall denote as ibi 
in this section, in different sizes of subvolumes. For the parameters we assume b\ — b <2 = 1 
and /nl = 50. We find that contributions from late-time evolution (black solid and dashed 
lines for *&z,,spt(&) and ib^fo (k), respectively) are redshift-independent, but the ones from 
primordial non-Gaussianity (colored lines at different redshifts) increase with increasing 
redshift. This is due to the redshift-dependence in M(k,a). Namely, while if>L,SPT(&) 
and ibL : b 2 (k) are independent of D(a), ib[ j j NL (k) is proportional to 1/ZD (a). This means 
that it is more promising to hunt for primordial non-Gaussianity in high-redshift galaxy 
surveys. We also find that for a given subvolume size ibLj NIj (k ) is fairly scale-independent, 
as 7 &l i spt(^) and ibL,b 2 {k). This is somewhat surprising because when k (the scale of the 
position-dependent power spectrum) is large we reach the squeezed limit, and this should 
be the ideal region to search for primordial non-Gaussianity. However, it turns out that 
what really determines the amplitude of ibj y j Nl (k) is the subvolumes size, as we can see from 
different panels in figure |2.1| One can understand this by considering the long-wavelength 
mode, ki , and the short-wavelength modes, k s (k s k s ). In the squeezed limit, 

B fNh (ki, k s , k a ) oc M(ki)M 2 (k s )[2P^(ki)Pq,(k s ) + P£(k s )] 
oc M(k l )M 2 (k s )P^(k l )P^(k s ) , 


(2.30) 
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as P<s>(k) oc k ns 4 with n s = 0.95. 
respectively, hence 

ih Lj^(k) oc 


Also M(k ) oc k 2 and k° on large and small scales 


n 3 —2 i,ns—4 


k i k 


s—2 




a 


oc 


a 


L,l 


(2.31) 


which is /^-independent and the amplitude of *&l,/ nl (^) is solely determined by ki , the 
subvolume size. This means that to hunt for primordial non-Gaussianity, it is necessary 
to use different sizes of subvolumes to break the degeneracy between ib L j^{k) and the 
late-time contributions. 


2.2 In configuration space 

2.2.1 Position-dependent correlation function 

We now turn to the position-dependent two-point statistics in configuration space, i.e. 
position-dependent correlation function. Consider a density fluctuation field, <5(r), in a 
survey (or simulation) volume V r . The global two-point function is defined as 


£0) = (h(x)h(x + r)) , 


(2.32) 


where we assume that <5(r) is statistically homogeneous and isotropic, so £(r) depends only 
on the separation r. As the ensemble average cannot be measured directly, we estimate 
the global two-point function as 


1 


d 2 1 




d 3 x h(r + x)<5(x) . 


'x.x+reV r 


The ensemble average of eq. (2.33) is not equal to £(r). Specifically, 


!«'» - E I £ 


' x,x+rGVr 


1 P d^v 

d 3 x (<5(r + x)<5(x)) =C(r)y / 


(2.33) 


d 3 x . (2.34) 


' x,x+r£l4 


The second integral in eq. (2.34) is V r only if r = 0, and the fact that it departs from V r is 
due to the finite boundary of V r . We shall quantify this boundary effect later in eq. (2.37). 

We now define the position-dependent correlation function of a cubic subvolume Vl 
centered at iq, to be 


((r,r L ) = - 


d 3 x 5(r + x)h(x) 


'x.r+xeVi 


= — y d 3 x A(r + x)h(x)W L (r + x - r L )W L (x - r L ) , 


(2.35) 


where Wl{ r) is the window function given in eq. (2.4). In this dissertation, we consider 
only the angle-averaged position-dependent correlation function (i.e. monopole) defined 
by 


d 2 f 


d 2 r 


€(r, r L ) = / — ^(r,r L ) = — — d 3 x h(r+x)h(x)W L (r+x-r L )W L (x-r L ) . (2.36) 


4(/r 


V l J 4tt 
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2. Position-dependent two-point statistics 


Again, while the overdensity in the entire volume 5 = J v d 3 r h(r) is zero by construction, 
the overdensity in the subvolume 5(r^) = J v d 3 r <5(r) = f v d 3 r h(r)lPi(r —r^) is in general 
non-zero. 


Similarly to that of the global two-point function, the ensemble average of eq. (2.36) is 
not equal to £(r). Specifically, 

A 1 f d 2 r f 

(f(r, r L )) = — I — I d 3 x (h(r + x)5(x))W L (r + x-r L )lF L (x-r L ) 

= £(r)^~ J I d?>x ' W l( v + x / )Vh L (x / ) = £(r)/ L , bndry (r) , (2.37) 

where ./nbndry(r) is the boundary effect due to the finite size of the subvolume. While 
/l, bndr y (' r ) = 1 for r — 0, the boundary effect becomes larger for larger separations. The 
boundary effect can be computed by the five-dimensional integral in eq. (2.37). Alterna¬ 


tively, it can be evaluated by the ratio of the number of the random particle pairs of a 
given separation in a finite volume to the expected random particle pairs in the shell with 
the same separation in an infinite volume. We have evaluated /l, tmdry(r) in both ways, and 
the results are in an excellent agreement. 

As the usual two-point function estimators based on pair counting (such as Landy- 


Szalay estimator which will be discussed in section 5.1.2) or grid counting (which will 


be discussed in appendix [C]) do not contain the boundary effect, when we compare the 
measurements to the model which is calculated based on eq. (2.36), we shall divide the 
model by /l, bndry(r) to correct for the boundary effect. 


2.2.2 Integrated three-point function 

The correlation between £(r, iq,) and 5(i\l) is given by 


(i(r, t l )S(t l )) 


V 2 I Itt" J d ’ ixi J (^( r + x i)£( x i)<5( x 2)) 

x W L (r + xi — r L )W L (xi - r L )W L (x 2 ~ r L ) 



x W L (r + xi)W L (xi)W L (x 2 ) , (2.38) 


where ('(r 1 ,r 2 ,r 3 ) = (5(r 1 )h(r 2 )5(r 3 )) is the three-point correlation function. Because of 
the assumption of homogeneity and isotropy, the three-point function depends only on the 
separations |r* — r ? for i ^ j, and so (£(r, ri)h(r^)) is independent of r^. Furthermore, as 
the right-hand-side of eq. (2.38) is an integral of the three-point function, we will refer to 
this quantity as the “integrated three-point function,” iQtX 7 ') = (£(r, r L )h(r L )). 

< L (r) can be computed if C(r!,r 2 ,r 3 ) is known. For example, SPT with the local bias 
model at the tree level in real space gives 


CM = &?CsPT(r) + blb 2 C b 2 (r) 


(2.39) 
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where £spt and fa 2 are given below. Here, b\ and 6 2 are the linear and quadratic (nonlinear) 
bias parameters, respectively. Because of the high dimensionality of the integral, we use 
the Monte Carlo integration routine in the GNU Scientific Library to numerically evaluate 
the eight-dimensional integral for ifaX 7 ')- 
The first term, Cspt, is given by [711 [12] 


Cspt(i~i, r 2 , r 3 ) 


-y6( r 12)6( r 23) + hl2,23[£«( r 12)0z(' r 23) + ( r 23)0Z ( r 12)] 

o ^( r 12 ) 0 K r 23 ) 6(n2)</>z( r 23) 66 ~ 23)^12) 



ri 2^13 


^2 3 


2 

^ 12,23 


6 (^ 12 ) + 


30K r 12) 


r 12 


6(r 23) + 


n 2 

30z( r 23) 

^2 3 


+ 2 cyclic , 


(2.40) 


where r i2 = | ri — r 2 |, /ii 2i23 is the cosine between ri 2 and r 23 , ' refers to the spatial derivative, 
and 

/ rlU f fjb 

— k 2 P l (k)sinc(kr), fa (r) = J — Pi(k)smc(kr ), (2.41) 

with Pi(k) being the linear matter power spectrum. The subscript l denotes the quantities 
in the linear regime. The second term, fa 21 is the nonlinear local bias three-point function. 
The nonlinear bias three-point function is then 


Cb 2 ( r i, r 2, r 3 ) = 6(^12)6^23) + 2 cyclic . (2.42) 

Note that £spt and fa 2 are simply Fourier transform of Hspt and B b2 respectively, as shown 
in eq. (2.26). 


Figure 2.2 shows the scale-dependencies of zCl.spt and i(,L,b 2 at z = 0 with Pi(k) com¬ 
puted by CLASS [STj. We normalize iQjfafa by erf;, where 


°L,l = (M*l) 2 ) = ^ 


d 3 k 


P,(*)|W^(k)| s 


(2.43) 


is the variance of the linear density field in the subvolume V^. The choice of this normal¬ 
ization is similar to that of the integrated bispectrum as discussed in section 2.1.2, and we 
shall discuss more details in section 2.2.4, We find that the scale-dependencies of ?Cl,spt(^) 
and i(L,b 2 ( r ) are similar especially on small scales. This is because the scale-dependence of 
the bispectrum in the squeezed limit is (see appendix of appendix |A.3[) 


-Bspt > 


68 ldlnk 3 Pi(k ) 
21 3 dink 


Pi(k)Pi(q ) , B b2 2 PXk)PXq) , 


(2.44) 


where k and q are the short- and long-wavelength modes, respectively. For a power-law 
power spectrum without features, the squeezed-limit £>spt and B b2 have exactly the same 
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2. Position-dependent two-point statistics 



Figure 2.2: Normalized i(l,spt (left) and i(i^b 2 (right) for L = 100 h 1 Mpc (red solid), 
200 h _1 Mpc (green dashed), and 300 h~ 3 Mpc (blue dotted) at z — 0. 


scale dependence and cannot be distinguished. This results in a significant residual degen¬ 
eracy between b\ and b 2l and will be discussed in chapter [5] where we measure the integrated 
three-point function of real data and fit to the models. When r is small, fCi( r )/ cr lz be¬ 
comes independent of the subvolume size. We derive this feature when we discuss the 


squeezed limit, where r L, in section 2.2.4 


2.2.3 Connection to the integrated bispectrum 

Fourier transforming the density fields, the integrated three-point function can be written 
as 


*C l(t ) = ^ 


d 3 q 1 


d 3 q 6 


V? J (2vr) 3 ,/ (2vr) 3 


(27r) 9 5 D (qi + q 2 + q 3 )Mqi + Q 2 + q4 + qs)fe(q3 + qe) 


X B{ qi , q 2 , q 3 )fFL(q4)hF L (q 5 )hFL(q6)e* [r -(qi+qd-r i -(q4+q5+q6)] 


d 3 k 

J2nf 


iB L (k)e 


irk 


(2.45) 


and it is simply the Fourier transform of the integrated bispectrum. Similarly, the angle- 
averaged integrated three-point function is related to the angle-averaged integrated bispec¬ 


trum (eq. (2.12)) as 


Kl[t) = j iB L (k)smc(kr) . 


(2.46) 


In chapter [5j we measure the integrated three-point function of real data, and thus we 
need the model for «Ci( r ) hi redshift space. Unlike the three-point function in real space 


(eq. (2.40) and eq. (2.42)), we do not have the analytical expression for redshift-space 
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Figure 2.3: Normalized ?Cl,spt (left) and ?Cl,6 2 (light) for L = 100 h~ l Mpc (red solid), 
200 h~ l Mpc (green dashed), and 300 h~ l Mpc (blue dotted) at z = 0 evaluated from 
Fourier space of a nine-dimensional integral (eq. (2.12) and eq. (2.46)). The same quantities 
evaluated from configuration space of an eight-dimensional integral (eq. ( |2.38| )) are shown 
with symbols. 


three-point function in configuration space. Since the integrated three-point function is 
the Fourier transform of the integrated bispectrum, we compute the redshift-space inte¬ 
grated three-point function by first evaluating the redshift-space angle-averaged integrated 


bispectrum with eq. (2.12 2 and then performing the one-dimensional integral eq. (2.46). 


This operation thus requires a nine-dimensional integral. 

To check the precision of numerical integration, we compare the results from the eight¬ 
dimensional integral in eq. (2.38) with the nine-dimensional integral eq. (2.12) and eq. (2.46) 
for both iCi,sPT and i(,L,b 2 - As the latter gives a noisy result, we apply a Savitzky-Golay 
filter (with window size 9 and polynomial order 4) six times, and the results are shown 
in figure 2.3. We find that, on the scales of interest (30 h^ 1 Mpc < r < 78 h^ 1 Mpc, 
which we will justify in section 5.1.3), both results are in agreement to within 2%. As 
the current uncertainty on the measured integrated correlation function presented in this 
dissertation is of order 10% (see section 5.2 for more details), the numerical integration 
yields sufficiently precise results. 


2.2.4 Squeezed limit 

In the squeezed limit, where the separation of the position-dependent correlation function 
is much smaller than the size of the subvolume (r <C L), the integrated three-point function 
has a straightforward physical interpretation, as for the integrated bispectrum. In this case, 

2 The explicit expression of the SPT redshift-space bispectrum is in appendix 


A.2 





















26 


2. Position-dependent two-point statistics 


the mean density in the subvolume acts effectively as a constant “background” density (see 
chapter [i] for more details). Consider the position-dependent correlation function, £(r, r^), 
is measured in a subvolume with overdensity h(r^). If the overdensity is small, we may 
Taylor expand £(r, r^,) in orders of 5 as 


€( r,r L ) 


£( r ) 1 , 5=0 + 


^(_r) 

dS 


6 + 0(P) . 

6 =0 


(2.47) 


The integrated three-point function in the squeezed limit is then, at leading order in the 
variance (h 2 ), given by 


<L(r) = (£(r,r L )5(r L )) = 


d£( r) 


dd 


(d 2 ) + 0(d 3 


(2.48) 


5=0 


As (d 2 ) = crfrl i<%(r) normalized by cr| is d£(r)/dd at leading order, which is the linear 


response of the correlation function to the overdensity. Note that in eq. (2.48) there is no 
dependence on the subvolume size apart from cr 2 , as shown also by the asymptotic behavior 


of the solid lines in figure 242 for r —» 0. 

As z<%(r) i s the Fourier transform of iB^k), the response of the correlation function, 
d£(r)/d5, is also the Fourier transform of the response of the power spectrum, dP(k)/dd. 
For example, we can calculate the response of the linear matter correlation function, 
d£i(r)/d5 , by Fourier transforming the response of the linear power spectrum, which is given 
in eqs. (2.23)—(2.24). In figure 2.4, we compare the normalized gCl.sptO") with d£i(r)/dd. 


Due to the large dynamic range of the correlation function, we divide all the predictions 
by £(r). As expected, the smaller the subvolume size, the smaller the r for ^Cl.sptM to be 
close to [1 /£i(r)][d£i(r)/dS], i.e., reaching the squeezed limit. Specifically, for 100 h Mpc, 
200 /r _1 Mpc, and 300 h^ 1 Mpc subvolumes, the squeezed limit is reached to 10% level at 
r ~ 10 hr 1 Mpc, 18 h~ x Mpc, and 25 h^ 1 Mpc, respectively. 


2.2.5 Shot noise 

If the density held is traced by discrete particles, A^(r), then the three-point function 
contains a shot noise contribution given by 


(<5d(ri)5d(r 2 )<5 d (r 3 )) = (<5(ri)<S(r 2 )5(r 3 )) 

'(h(ri)h(r 2 )) 


n(r 3 ) 


M r i - r 3 ) + 2 cyclic 


+ 


M r ! - r 2 )<dD(ri - r 3 ) 


n(r 2 )n(r 3 ) 


(2.49) 


where n(r) is the mean number density of the discrete particles. The shot noise can be 
safely neglected for the three-point function because it only contributes when ri = r 2 , 

3 If 6 = 5i then a 2 L = o\ l . But 5 can in principle be nonlinear or the mean overdensity of the biased 
tracers, so here we denote the variance to be 
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Figure 2.4: The linear response function [1 /£i(r)][d£i(r)/dS] (black solid) and the nor¬ 
malized *CL,SPT( r ) for L = 100 h~ l Mpc (red dotted), 200 h~ l Mpc (green dashed), and 
300 h l Mpc (blue dot-dashed). The light and dark bands correspond to ±5% and ±10% 
of the predictions, respectively. 


r 3 = r 3 , or r 2 = r 3 . On the other hand, the shot noise of the integrated three-point 


function can be computed by inserting eq. (2.49) into eq. (2.38), which yields 


<shot(r) =£( r )^2 / ^ / d3 ' x 


+ 


_n(x±r±r L ) n(x ± r L ) 


W L (x + r)Wi(x) , (2.50) 


where we have assumed r > 0. If we further assume that the mean number density is 
constant, then the shot noise of the integrated three-point function can be simplified as 


iC.hc.tM = 2CW— /i.bndryM • 


(2.51) 


For the measurements of PTHalos mock catalogs and the BOSS DR10 CMASS sample 
shown in chapter |5j the shot noise is subdominant (less than 7% of the total signal on the 
scales of interest). 
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Chapter 3 

Separate universe picture 


Mode coupling plays a fundamental role in cosmology. Even if the initial density fluctua¬ 
tions generated by inflation are perfectly Gaussian, the subsequent nonlinear gravitational 
evolution couples long and short wavelength modes as well as generates non-zero bispec¬ 
trum. A precise understanding of how a long-wavelength density affects the small-scale 
structure formation gravitationally is necessary, especially for extracting the signal of bis¬ 
pectrum due to primordial non-Gaussianity. 

A useful way to describe the behavior of the small-scale structure formation in an 
overdense (underdense) environment is the separate universe picture [ITT], 11121 150 i :TTf 11471 
EH ESI Sal EH- Imagine a local observer who can only access to the comoving distance 
of the short-wavelength modes ~ 1 /ks- If there exists a long-wavelength mode such that 
1/fei > 1 /ks, then the local small-scale physics would be interpreted with a Friedmann- 
Lemaitre-Robertson-Walker (FLRW) background [93]. That is, in the separate universe 
picture, the overdensity is absorbed into the modified background cosmology, and the 
small-scale structure formation evolves in this modified cosmology. 

The separate universe picture can be proven from the general relativistic approach: one 
locally constructs a frame, Conformal Fermi Coordinates, such that it is valid across the 
scale of a coarse-grained universe A -1 (Aq, < A < ks) at all times [12311431 144] , In Conformal 
Fermi Coordinates, the small-scale structure around an observer is interpreted as evolving 
in FLRW universe modified by the long-wavelength overdensity. The separate universe 
picture is restricted to scales larger than the sound horizons of all fluid components, where 
all fluid components are comoving. 

In the presence of 5, the background cosmology changes, and the position-dependent 
power spectrum is affected accordingly. Since the response of the position-dependent power 
spectrum to the long-wavelength overdensity, [dP(k)/dS] t 5 =0 , can be related to the bispec¬ 
trum in the squeezed configurations, the separate universe picture is useful for modeling 
the squeezed-limit bispectrum generated by nonlinear gravitational evolution. Not only 
this is conceptually straightforward to interpret, but it also captures more nonlinear effects 
than the direct bispectrum modeling, as we will show in chapter [4j 

In this chapter, which serves as the basis for the later separate universe approach mod¬ 
eling, we derive the mapping between the overdense universe in a spatially flat background 
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3. Separate universe picture 


cosmology (but with cosmological constant) to the modified cosmology in section 3.1 In 
section 3.2, we show that if the background cosmology is Einstein-de Sitter, i.e. matter 


dominated, then the changes in the scale factor and the linear growth factor can be solved 
analytically. 


3.1 Mapping the overdense universe to the modified 
cosmology 

Consider the universe with mean density p(t) and a region with overdensity 5 p (t), then the 
mean density p(t) in this region is 

p(t) = p(t)[l + 5 p {t)\ . (3.1) 

In this section, we shall derive the mapping of the cosmological parameters between the 
fiducial (overdense) and modified cosmologies as a function of the linearly extrapolated 
(Lagrangian) present-day overdensity 

Sa = s ^m’ (3 ' 2) 

where D is the linear growth function of the fiducial cosmology, to is the present time, and 
tj is some initial time at which 5 p (ti) is still small. However, we shall not assume Slo to be 

small. 

The mean overdensity of the overdense region can be expressed in terms of the standard 
cosmological parameters, i.e. p(t 0 ) = and H 0 = 100 h km s -1 Mpc -1 , as 

n m h\ c , .. & m h 2 , \ 

M l+m] = W) (33) 


where the parameters in the modified cosmology are denoted with tilde, 
cosmology, we adopt the standard convention for the scale factor such 
In contrast, for the modified cosmology, it is convenient to choose a(t 
S p (a —>■ 0) = 0, which then leads to 

For the fiducial 
that a(to) = 1. 
—» 0) = a and 

f 'Imh 2 = fl m h 2 . 

(3.4) 

We also introduce 


a(t) = a(t)[ 1 + 6 a (t)] , 

(3.5) 

so that 


1 + = [1 + ^a(t)] 3 • 

(3.6) 
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Using eq. (3.5), the first and second time derivatives, represented with dots, of the scale 
factor are 


a _ _ a[l + *5q] + aS a _ H + S a 

a a[l + 5J 1 + 5 a 

a _ a[l + 5 a \ + 2aS a + a5 a _ a | 5 a + 2HS a ^ 

__ r 1 I -It'* / 

a gl\\ + o a J a 1 -f- d a 

The two Friedmann equations in the flat fiducial cosmology are 

(£)’ =H\t) = & f-m+px(t)], 

- =~m + px(t) + 3pjt)], (3.8) 

CL O 

where px and p x are the energy density and pressure for dark energy, respectively. For the 
modified cosmology, the Friedmann equations hold, but with non-zero curvature K and 
modified densities and scale factor as 

• 2 ~ 

(a\ ~n, . 87rG r ~. . „ . K 

(aj = + 

CL 4:71 (sr v ~ x , \ ^ . ... , x 

V =-+ px(t) + 3 p x (t)] . (3.9) 

Before we derive the relation for the curvature K , let us first discuss the dark energy 
component. 

If dark energy is not a cosmological constant, then there are also perturbations in dark 
energy fluid, i.e. Spx = Px ~ Px and Spx = px ~ Px- In order for the separate universe 
approach to work, matter and dark energy have to be comoving and follow geodesics of 
the FLRW metric. Since this requires negligible pressure gradients, the separate universe 
approach is only applicable to density perturbations with wavelength 2ir/k that are much 
larger than the dark energy sound horizon, k <C H 0 /\c s \, where the sound speed is defined 
as (? s = fipx/bpx [$>]• This means that the region with the overdensity that we consider 
here has to be much larger than the dark energy sound horizon. For simplicity, in the 
following we shall assume that dark energy is just a cosmological constant A, thereby 
Pa = Pa and Pa = Pa = ~Pa- 

In order to be a valid Friedmann model, the curvature has to be conserved, i.e. K = 0. 
To show this, we express the curvature using the first Friedmann equation and take the 
time derivative: 


167T G_x.r~ „ , 87 tGU 2 ~ 167rG_^ r ~ „ . 247tG^^~ 

K = —-— aa[p + px\ -\ - a p — 2aa = - aa[p + paJ- o,ap — 2aa 


3 

87 tG 


aa\—p + 2 oa 1 — 2 aa = 2 aa- — 2 aa = 0 , 

a 


3 


(3.10) 
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where we use = 0 and the continuity equation p = —3pH. 

To solve the curvature of the modified cosmology in terms of the fiducial cosmological 
parameters and Slo, we express K by the difference of the first Friedmann equation between 
the modified and fiducial cosmologies as 


K 


87 tG 


p5 p {l + 5 a y — 2HS a (l + 5 a) — S t 


3 pv 

——P [(1 + 5a) 1 — (1 + ^a) 2 ] — 2HS a {l + 8 a ) — 8l , 


(3.11) 


where we use eq. (3.6) and eq. (3.7) to represent S p and H, respectively. Since the cur¬ 
vature is conserved (eq. (3.10)), eq. ( |3.11 ) can be evaluated at an early time £* at which 
the perturbations S a (ti) and S p (tj) are infinitesimal as well as the universe is in matter 
domination. In this regime, we have 


H 2 = Hln m a~ 3 


5 a = —5p/3 , 5 a = HS a , 


(3.12) 


with which we can derive 


k 

a 2 (ti ) 


87 tG 

~Y~ 


■p(ti)[-38a(ti)] 


2 H 2 (U)5 a (U) 


5 n m H^6 p (ti) 
3 a 3 (U) 


(3.13) 


This then to express K in terms of the fiducial cosmological parameters and Slo as 


K 5 Q, 


HI 


3 a{t/) 


Sp{U) = o 


5 a 


3 D{to 


JLO 


(3.14) 


where we use the fact that in the matter dominated regime D{t/) = a{t/). 

We now derive the cosmological parameters, Q m , 12 a, and Cl K of the modified cosmology. 
Note that by convention these parameters are defined through the first Friedmann equation 
at to where a(t 0 ) = 1 and H(to) = H 0 , therefore 


n K = — A 

Ft 2 

n o 


~ 87 tG ~ ~ 

“m = P{to) 


~ 87 tG 

nA = Wf A 


(3.15) 


The cosmological parameters in the modified cosmology can be related to the ones in the 
fiducial cosmology through H 0 = H 0 (l + Sh), which become 


&m — ^m(l + <5ff) 2 , — ^a(1 + £u) 2 , — 1 7l m Ha — 1 (lT^/f) 2 , (3.16) 


where we use Vt m + Ha = 1 because the fiducial cosmology is flat. 
Alternatively, Sr can be expressed in terms of K as 



- 1 . 


(3.17) 
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There is no solution for 5h if K/H q > 1, or equivalently Slo > (|d(^)) _ 1 - This is because 
for such a large positive curvature, the universe reaches turn around before a — 1, at which 
the modified cosmological parameters are defined. In other words, this is not a physical 
problem, but merely a parameterization issue under the standard convention. 

Finally, we shall derive the equation for 5 a (t), so that the observables of different cos¬ 
mologies can be compared at the same physical time t. Inserting the second Friedmann 


equation of the fiducial and modified cosmologies into eq. (3.7) yields an ordinary differ¬ 
ential equation for the perturbation to the scale factor 

4v tG 


5 a + 2 H5 a H —-—p [(1 + 6 a ) 2 — (1 + <5 a )] — 0 , 


or equivalently 


, , 4 S 2 0 

5 0 + 2 H5 a — —- 


3 1 + 8, 


- AirGp5 p (l + S p ) — 0 . 


(3.18) 


(3.19) 


When linearizing (8 P -C 1) eq. (3.19), one obtains the equation for the linear growth factor. 
More generally, eq. ( |3.19[ ) is exactly the equation for the interior density of a spherical 
top-hat perturbation in a ACDM universe [138] . For a certain f out , one can numerically 


calculate <5 a (f 0 ut) through eq. (3.18) to get a(f out ) = a(f out )[l + <5 a (f ou t)]- Alternatively, one 
can numerically evaluate a(f ou t) by 


'“ttpout) 


taut. 


da 


'‘A (t-OUt) 


da 


ra(tout ) [14“<5a (tout)] 


da 


aH(a) 


aH(a) 


aH(a) 


(3.20) 


3.2 


The modified cosmology in Einstein-de Sitter back¬ 
ground 


In the Einstein-de Sitter (EdS) universe we have 

/3 \ 2//3 

= 1 ; Pa = 0 ; H{a) = H 0 a~ 3/2 ; a(t) = ( -H 0 t j ; D(t) = a(t) . (3.21) 

In this section, we shall show that if the background cosmology is EdS, then the scale 
factor S a (t) and the Eulerian overdensity S p (t), as well as the linear growth factor D(t) 
and the logarithmic growth rate /(f) in the modified cosmology (overdense region) can be 
expressed in series of ^o- 


3.2.1 Scale factor and Eulerian overdensity 

In order to solve S a (t), we consider the spherical collapse model for a overdense region 
[66UI2SUI22I. Since the scale of this region is much smaller than the horizon size, we can 
use the Newtonian dynamics and write the equation of motion of a shell of particles at 
radius r as 

» GM ~ Ait n . 
r = -—— , M = — p{ti)r (ti) = constant , 


(3.22) 
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where U is some initial time. Note that in eq. (3.22) we neglect the shell crossing, i.e. if 
fi(tj) > f 2 (ti) then fi[t) > f 2 (t) for all t. 


Eq. (3.22) is known as the “cycloid” and the parametric solution is 

d(9) = —— = A{1 — cos 9) , t{9) = B{9 — sin 6) , 

x 


(3.23) 


where x is some comoving distance for the normalization. Using the Leibniz rule, one finds 
that 

dr d6 Ax sin 6 ■■ Ax 1 

1-COS0) ' r = -“- (3 ' 24) 

and the equation of motion thus requires 


B 2 (1 — cos 6) 2 


7L 3 x 3 ~ AnG ,_ 3 , ._ 3 A 3 

_ =GM = —p( U )a 3 (U)x 3 ^ 


(3.25) 


Note that since this region is overdense, Q m — 1 = —Qk > 0 and it is positively curved. 
For an underdense (negatively curved) region, the similar parameterization works with cos 
and sin in eq. (3.24) replaced by cosh and sinh. 

To determine the constants A and B, we can write the first Friedmann equation using 
the parametric solution as 

2 i n _ _2 m 

1-3 


TTZ = I Z , = 1 (1 COS 2 9) _ ^2 

a ) B 2 (1 — cos#) 4 


(1 - cos 2 9) _ Hi 
B 2 “ A 2 


^-(sL-i) 


= Hi \ fl m [A(l - cos #)] + (1 - O m )[A(l - cos i 

a 


1-2 


— cos 9 


-j- - 2(n m -1) 


— cos #(O m — 1) 


(3.26) 


Since eq. (|3.26|) is valid for all 6, the coefficient multiplied by cos 9 must be zero, which 

11 in 


then allows us to solve 


A- 1 - 


2 (Sin ~ 1) 


B = i- 


a 


2 H 0 (Q m — l) 3 / 2 


(3.27) 


One finds that eq. (3.27) indeed yields A 3 /B 2 = H'lVt m /2. 

For simplicity, we define 

^ m ~ 1 1 /-i i x \ 2 ^ ^ x 

e - —p.-- 1 - (1 + o H ) - 772 “ , 

i L rrii n o 


so that 


a{9) — X (1 — cos 9) ; t{9) = H 0 e 3//2 (# — sin#) 


(3.28) 

(3.29) 


The goal is to obtain 

a(t) = a(£)[l + S a (t)] = a(t) 


1 + y^e w [a(t)(5x,o] 


n= 1 


a (0) — 1 + ^2 e n^L0 ) (3.30) 


n=1 
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we thus need to solve 0 q such that 


*(0o) =t 0 = ^t 0 e 3/2 (0 O - sin 6*o) B 0 - sin0 o = 3/2 

We perform a series expansion, 


(3.31) 


11 

0 O - sin 0 O = ^0o - J7^ e t + ' ’ ’ = bn6 ' 

n =1 


2n+l 

0 


6 u 120 

and solve 0 q order by order. That is, for the n th order solution, we solve 


(3.32) 


i e 3 /2 = 




(n— fc+1) 


1 2fc+l 


k= 1 


(3.33) 


and 0 O = lim^oo 0q‘\ Note that in order to trust the final solution at order <5™ 0 , the series 
solution needs to be expanded to n — m + 1. In the following we choose m — 5, which 
yields 

1 2 9 4 o 43 

- e 2 -A - e 6 - 4 - - 

175 


0n = 2c 1 / 2 


1 + — e + —e 2 + —^-e 3 + —^-e 4 H- 


15 


1575 


67375 


(3.34) 


Finally, we insert eq. (3.34) into 5(0), expand in e, replace e with an d match with 


eq. (3.30) for the coefficients e n . For the first five coefficients, we get 

1894 


er - -3 ; e 2 - ; e 3 - 


23 

1701 


3293 


; e4 — 


392931 


; 65 


1702701 


(3.35) 


Once S a (t) is known, we can use eq. (3.6) to solve the Eulerian overdensity 5 p (t) = 
fn[3TA)a(t)] n , and the first five coefficients are 


, _ f _ 17 . _ 341 . _ 55805 _ 213662 

Ji- , h - 21 , h - 5g7 , J4 - 130997 > h ~ 729729 


(3.36) 


While eqs. (3.35)—(3.36) are strictly valid only if the background is EdS, we find that 
they are also accurate in ACDM background if a(t) is replaced with D(t)/D(to), where 
D(t) is the growth factor in the fiducial cosmology. In other words, 


OO 

71—1 


m = E 


L 0 


m 

D(t„) 




71—1 


L0 


m 

D(t 0 ) 


(3.37) 


Figure 3.1 quantifies the performance of the EdS expansion with a(t) being replaced by 
D(t)/D(t 0 ) for S a (t), i.e. eq. (3.37). For a large range of 5lo, the results of EdS expansion 
agree very well with the numerical solution of the differential equation eq. (3.18). Specif¬ 
ically, at the third (fifth) order expansion, the fractional difference is at the sub-percent 
level for |<5z, 0 | ~ 1- 




























36 


3. Separate universe picture 



1.10 


1.05 


O 

o 

C 

O 

'in 

C 

fD 

Q. 

X 

<D 


1.00 


0.95 


0.90. 


- 1.0 


expansion 1st 
expansion 3rd 
expansion 4th 
expansion 5th 



Figure 3.1: (Left) Perturbation in the scale factor in the modified cosmology, S a (t 0 ), as 
a function of 8lo- The present-day scale factor is a(t 0 ) = 1 in the background ACDM 
universe, with Vt m = 0.27 and 12 a = 0.73. The black solid line shows the numerical solution 


to the ordinary differential equation, eq. (3.18); the yellow solid, green dot-dashed, blue 


dotted, and red dashed lines show the series solutions in EdS, eq. (3.37), at the first, third, 


fourth, and fifth order, respectively. (Right) Same as the left panel, but for the ratios of 
the series solutions expansion to the solution of the ordinary differential equation. 


3.2.2 Small-scale growth 

In this subsection, we shall iteratively so) 
growth rate f(t ) in series of 5l 0 , as eqs. (|3.35[)-(|3.36[), in the modified cosmology with the 


ve th e linea r growth factor D(t) and logarithmic 


EdS background. We consider a small-scale (short-wavelength) density perturbation in the 
overdense region such that 5 S = p/p — 1. Note that S s should not be confused with the 
long-wavelength density perturbation S p of this entirely overdense region with respect to 
the EdS background. Moreover, 5 S is defined with respect to the background density of 
the overdense universe p instead of the EdS background p. 

The small-scale growth equation for S s in the modified cosmology is given by 


S s + 2 H8 S - 4t xGp5 s = 0 . 


(3.38) 


Using 


H 2 

AnGp 


Hi 


- - n i (1 


-3 


-O H 2 n ~ 3 
2 1 


+ (1 — Q m )a 




) 


(3.39) 
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the growth equation can be rewritten as 


S s + 2H 0 





0 . 


(3.40) 


Note that although we map the overdense universe to a positively curved universe, the 
curvature contribution to the Poisson equation is neglected. If K/ ~ 1, the correction 

in the Poisson equation becomes relevant for the small-scale modes 5 S that are around the 
horizon size. Since we are studying the subhorizon evolution of the small-scale modes, 
and moreover we are mostly interested in the overdensity such that K/H q ~ 5lq <C 1, 
the correction is entirely negligible. Thus, the curvature contributes to the growth only 
through the expansion rate H. 

Replacing the time coordinate t with y = lna(t) where a(t ) is the scale factor in the 
EdS background, we rewrite the growth equation as 


d 2 ? 

W‘‘ + 


/ r. \ 1/2 

2(1 + 5 a ) 3//2 (1 — -<5lo[ 1 + <5a] J — 


l + S a )- 3 l = 0 , (3.41) 

dy 2 


Thus far, all the derivations are exact. To see that eq. (3.41) makes sense, we consider the 
zeroth order approximation, i.e. > 0. In this regime, 5 a —* 0 and the growth equation 
becomes 

d^^ 0) + 2 dy^ ~ 2 ^ 0) = ° ’ ^ 3 ' 42 ) 

where the subscript (0) denotes that it is the zeroth order solution. There are two solutions 
to 5s , the growing mode proportional to a and the decaying mode proportional to a -3 / 2 . 
As expected, because 5lo —> 0, the result is identical to the growth in the background 
EdS cosmology. In the following, we shall drop the decaying mode following the standard 
practice. Furthermore, we shall normalize 5s to a(t) at early times, and replace it with 
D(t) to denote the small-scale growth factor. This means D^(t) = a(t). 

To solve higher order solutions, we insert the expansion of 5 a in terms of a{t)5Lo = e y 5Lo 


(eq. (3.30) and eq. (3.35)) into the growth equation and obtain 


d 2 

w D{y) + 


E 

171=0 


rm my 




E 

m =0 


A Xm my 
Uj m u 


D(y) = 0 , 


(3.43) 


with coefficients c m and d m given by 


/ 5 \ 1 / 2 s °° 

2(1 + 5 a ) 3 ^ (1 — -<5lo[ 1 + <^a] j — g = c m [a(t)5Lo\' 

^ ' m= 0 

O oo 

-(1 + <5a)~ 3 =^2d m [a(t)5 L o] 
m =0 


(3.44) 
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Correspondingly, we write the pure growing-mode solution in series of 8lq as 


D = J2 9 n5 


n p(n+l)y _ 
L0 e 




(3.45) 


n =0 


n =0 


with coefficients g n . Given our normalization, i.e. = a(f), go = 1. Thus, 

i oo oo 

D(y) = $> + = 5> + l) 2 g„^e‘” +1 >» 

y n =0 y n =0 


(3.46) 


Supposed that we have solutions of D(y) to the (n — l) th order, then the solution at 
the n th order has to satisfy 

n 

(n + 1 ) 2 gJ%e^ y + ^ g n _ m 8 n L - m e^~ m+ ^ y [(n - m + l)c m - d m ] = 0 . (3.47) 

m =0 

The term S^ 0 e^ n+1 ^ y factors out, and we obtain a simple algebraic relation for g n in terms 
of c m and d m for 0 < m < n, and g m for 0 < m < n — 1 as 


(n + 1 ) 2 g n + ^2 9n-m[(n - m + 1 )c m - d m \ = 0 . 


(3.48) 


m =0 


Using e„ in eq. (3.35) to get c m and d m through eq. (3.44), we obtain the first five g n to 
be 


9 1 = 


13 

21 


g 2 = 


71 

189 


9 3 = 


29609 

130977 


9 4 = 


691858 

5108103 


fib = 


8682241 

107270163 


(3.49) 


Similar to the previous subsection, the expansion of D in terms of 5lo with the coeffi¬ 


cients (eq. (3.49)) is strictly valid in the EdS background cosmology. In order to generalize 


from EdS to other cosmologies, we replace a(t) with D(t)/D(t 0 ), so that 


m = D(t) 


OO 


+ 2^,9* 


77.—1 


JLO 



(3.50) 


Figure [3)2| quanti fies th e performance of the expansion with a(t) being replaced by D(t)/D(t 0 ) 
for D(t), i.e. eq. ( |3.50 ). The agreement is not as good as for 5 a (t), nevertheless the fifth or¬ 
der expansion gives 5% fractional differences at \Slo\ ~ 1. Note also that while for positive 
6lo the EdS expansion is always smaller than (but converging to) the numerical solution 
in ACDM, for negative 5^0 the fourth order expansion has a different trend compared to 
the other orders. This is because at the n th order <5£ 0 dominates when \Slo\ > U an d at 
the negative 5^0 en d the result would thus depend on the parity of the expansion order. It 
is clear though the higher the expansion order, the better the agreement. 
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Figure 3.2: (Left) The linear growth factor in the modified cosmology, D(t 0 ) with a(t 0 ) = 1, 
in the background ACDM universe with Vt m = 0.27 and Da = 0.73. The black solid line 


shows the numerical solution to eq. (3.40); the yellow solid, green dot-dashed, blue dotted, 


and red dashed lines show the series solutions, eq. (3.50), at the first, third, fourth, and 


fifth order, respectively. (Right) Same as the left panel, but for the ratios of the series 
solutions to the solution of the differential equation of the small-scale growth. 


With the expansions of 8 a (t) (eq. (3.37)) and D(t ) (eq. (3.50)), we can finally derive 
the series expansion for the logarithmic growth rate, 


/(*) = 


d\nD(t ) D(t) a(t) 


d\na(t) D(t) a(t) 

In the modified cosmology, we have (defining e$ = 1 as for g 0 ) 

7, t s = Da = DT,n=^ n + l )9nK( t ) _ a Y.n =o eJljt) _ 

Dd D Y,n=i)dJl{t) a^“ =0 e„5£(f)+ ag^“ =0 ne n (5£(f) 

= / E^o(n + i )g n m 


(3.51) 


E“=0 9nSl(t) J2n =0 e n 61(t) + f E“=o neJlit) ’ 


(3.52) 


where 8Lit) = 8lo~§^- Note that in the EdS background / = 1, and so eq. (3.52) can be 
simplified as 

7 (f] = So (n + 1 )gJljt) E,T=o e ^lit) ( „ 

1 ZZo 9nSl{t) Er=o(« + l)^W ' j 

Figure 3.3 shows the performance of the expansion of / in the ACDM background. It is 
not as good as for D, but for |5x, 0 | ~ 0.8 the fifth order expansion gives approximately 
5% fractional difference with respect to the numerical solution of the differential equation. 
One also notes the parity-feature at the negative 8lo, as for D. 
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Figure 3.3: (Left) The logarithmic growth rate in the modified cosmology, /(to) with 
a(to) = 1, in the background ACDM universe with Q m = 0.27 and = 0.73. The black 
solid line shows the numerical solution to eq. (3.40); the yellow solid, green dot-dashed, 
blue dotted, and red dashed lines show the series solutions, eq. (3.52), at the first, third, 
fourth, and fifth order, respectively. (Right) Same as the left panel, but for the ratios of 
the series solutions to the solution of the differential equation of the small-scale growth. 











Chapter 4 

Measurement of position-dependent 
power spectrum 


As introduced in section 2.1 
trum, 


the correlation between the position-dependent power spec- 


p(k, r J = y l / * (k “ qi),5 (“ k - ,q,+q,) , 

and the mean overdensity, 

* (r ‘> = V L / fly , 


is the integrated bispectrum, 

= (P(/c,r L )h(r L )) 

1 I" cPk I" d 3 qi 

V‘l J 4:71 J (27r) 3 


d 3 g 3 

(2vr) 3 


■ B( k - qi, —k + qi + q 3 , -q 3 ) 
xW^(q 1 )W i (-q 1 -.q 3 )W L (q 3 ) 


1 


(4.1) 


(4.2) 


(4.3) 


with Vl being the size of the subvolume. 

In the squeezed limit where the scale of the position-dependent power spectrum is much 
smaller than the subvolume size, i.e. ^ > 1/L, the integrated bispectrum can be simplified 
as 

iB L (k) “5” A / (^ W|(q )P(q)f[k)P(k) = olf(k)P(k) . (4,4) 

where f(k) = 2/(0, k) with /(ki,k 2 ) being a dimensionless symmetric function for the 
separable bispectrum, and a\ is the variance of the density fluctuation in Vl, 

a ‘ L = viS (dy ■ 


(4.5) 
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4. Measurement of position-dependent power spectrum 


An intuitive way to arrive at eq. (|4.4|) is to consider the expansion of the position-dependent 

(4.6) 


power spectrum in the presence of a long-wavelength density fluctuation 5 as 

dP(k ) 


P(k,i L ) = -P(*0| J=o + 


dl 


5 + • • • , 


5=0 


and the leading-order correlation between P(/c,r 0 ) and 5 is 


iB L (k) = 


2 d In P{k) 


dd 


p(k). 


(4.7) 


5=0 


Inspired by eq. (4.4) and eq. (4.7), we define the normalized integrated bispectrum to be 

, (4.8) 

P(k)al 

and it is equal to the linear response function, f{k) or din P(k)/dd, in the limit of kL —>- oo. 

In this chapter, we measure the position-dependent power spectrum and the integrated 
bispectrum from IV-body simulations in section |4.1 and compare with the theoretical 
modeling of the measurements in section [42] and section 43 At the end of this chapter, we 
shall discuss the dependence of the integrated bispectrum on the cosmological parameters 
in section 4.4, and the expected constraint on the primordial non-Gaussianity using the 
Fisher matrix calculation in section 14751 We conclude in section 14761 


4.1 iV-body simulations and the estimators 

We now present measurements of the position-dependent power spectrum from 160 colli¬ 
sionless Wbody simulations of a 2400 h~ l Mpc box with 768 3 particles (which corresponds 
to 2.29 x 1O 12 M 0 ). The same simulations are used in [47], and we refer to section 3 of [47] for 
more details. In short,the initial conditions are set up using different realizations of Gaus¬ 
sian random fields with the linear power spectrum computed by CAMB [96], [95]. We adopt 
a flat ACDM cosmology, and the cosmological parameters are Q m = 0.27, f \h 2 = 0.023, 
h = 0.7, n s = 0.95, and cr 8 = 0.7913. The particles are displaced from the initial grid 
points using the second-order Lagrangian perturbation theory [40] at the initial redshift 
Zi = 19. The simulations are carried out using the Tree-PM code Gadget-2 Dsa. taking 
only the gravitational force into account. 

To construct the density fluctuation field on grid points, we first distribute all the 
particles in the 2400 h Mpc box onto a 1000 3 grid by the cloud-in-cell (CIC) density 
assignment scheme. Then the density fluctuation held at the grid point r g is 

i(r„) = _ ! , (4.9) 

where hat denotes the estimated quantities, N(r g ) is the fractional number of particles 
after the CIC assignment at r g , and N = 768 3 /1000 3 is the mean number of particles in 
each grid cell. 
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We then divide the 2400 h -1 Mpc box in each dimension by N cut = 4, 8, and 20, so that 
there are 64, 512, and 8000 subvolumes with a side length of 600, 300, and 120 h -1 Mpc, 
respectively. The mean density perturbation in a subvolume centered at is 

*( r i) = ^ T a) ’ ( 4 - lQ ) 

V id rg ev L 


where (lV gr id) 3 = (1000/lV cut ) 3 is the number of grid points within the subvolume. To 
compute the position-dependent power spectrum, we use FFTfj^jto Fourier transform <5(r ff ) 
in each subvolume with the grid size (lV gr id) 3 . While the fundamental frequency of the 
subvolume, k F = 2i t/L, decreases with the subvolume size L, the Nyquist frequency of the 
FFT grid, k Ny = k F N gvid /2 ~ 1.3 h Mpc -1 , is the same in all cases. 

The position-dependent power spectrum is then computed as 


P(k, r L ) 


1 

Vf, N mode 


i^( k i> r i)i 2 > 

k-Ak/2<\'k i \<k+Ak/2 


(4.11) 


where lV mo de is the number of Fourier modes in the bin [k — Ak/2, k + Ak/2], and we set 
A k ~ 0.01 h Mpc -1 in all cases. We choose this Ak for all N cut to sample well the baryon 
acoustic oscillations (BAO) and thereby are able to show how the window function of the 
different subvolumes damps the BAO (see figure 4.2). We follow the procedures in 


to correct for the smoothing due to the CIO density assignment and also for the aliasing 
effect in the power spectrum. Note, however, that this correction is only important for 
wavenumbers near the Nyquist frequency 1.31 h Mpc -1 , and we are interested in scales 
k < 0.4 h Mpc -1 . 


Figure |4.1| shows the position-dependent power spectrum measured from 512 subvol¬ 
umes with L = 300 h -1 Mpc in one realization at z = 0. The color represents S(r F ) of each 

subvolume. The positive correlation between the subvolume power spectra and S(r L ) is 
obvious. The response of the position-dependent power spectrum to the long-wavelength 
density fluctuation is clearly measurable at high significance in the simulations. 

We measure the integrated bispectrum through 


NT 3 

iV cut 


iB L {k ) = 


Nl 




P(k,r Lti )S(r Lti 


(4.12) 


cut 


i— 1 


where P(k, rp ti ) and S(r F i ) are measured in the i th subvolume. Further, motivated by 
eq. (4.4), we normalize the integrated bispectrum by the mean power spectrum in the 
subvolumes, 


^c 3 ut 


* ’cut 


(4.13) 


2=1 


■^ast Fourier Transformation library: www.fftw.org 
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Figure 4.1: Position-dependent power spectra measured from 512 subvolumes with L = 
300 h Mpc in one realization at z = 0. The color represents 8(tl) of each subvolume. 


and the variance of the mean density fluctuation in the subvolumes, 

TV 3 
1 v 


*L = 


i 


N* 




(4.14) 


cut •_ 


i=l 


Note that by construction 


This quantity 


nL 


8l = 


N 3 + 

1 'cut 




= 0 


i— 1 


iB L (k) 

P L (k)&l 


(4.15) 


(4.16) 


is the estimator of the normalized integrated bispectrum (eq. ( |4.8 )), and is equal to the 
linear response function, dln.P(k)/dS, given in eq. (4.7) in the limit of kL —y oo. 


Figure 4.2 shows the normalized integrated bispectrum, averaged over 160 collisionless 


IV-body simulations at different redshifts. For clarity, no error bars are shown in this figure. 
We have compared the results with a higher-resolution simulation with 1536 3 particles and 
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Figure 4.2: Normalized integrated bispectrum averaged over 160 collisionless IV-body simu¬ 
lations with Gaussian initial conditions. From left to right are N cut =4 (L — 600 h~ l Mpc), 
8 (300 h^ 1 Mpc), and 20 (120 h^ 1 Mpc); the blue, green, yellow, and red lines are z — 3, 
2, 1, and 0, respectively. For clarity, we do not show the error bars. 


starting at higher redshift (z{ = 49 compared to z t = 19 for our 160 simulations). For 
the scales and redshifts shown in figure 42, the differences are less than 1%. However, we 
expect an up to 5% uncertainty in the integrated bispectrum at z = 3 (less at lower z) due 
to transients which affect the bispectrum more strongly than the power spectrum |40 lllll| . 
as well as other systematics such as mass resolution. 

Since the initial conditions are Gaussian, the bispectrum is generated entirely by non¬ 
linear gravitational evolution. We thus measure the effect of a long-wavelength density 
perturbation on the evolution of small-scale structures. The wiggles visible in each panel 
of figure 42 are due to the BAOs. The BAOs in the right panel are strongly damped 
because the box size (120 h Mpc) approaches the BAO scale, and the window function 
smears the BAO feature |30j . Further, BAO amplitudes are larger at higher redshifts as 
they are less damped by nonlinear evolution [53]. The broad-band shape of the normalized 


integrated bispectrum evolves on small scales due to nonlinear evolution, leading to an 
effective steepening of its slope. We now turn to the theoretical modeling of the results 


shown in figure 4.2 


4.2 Bispectrum modeling 

We use two different approaches to model the integrated bispectrum. In the first approach, 
we model the bispectrum and compute the integral to obtain the integrated bispectrum. 
In the second approach, we model the response of the small-scale power spectrum to a 
long wavelength perturbation directly using the “separate universe” picture. For clarity, 
we will show the comparison between model prediction and simulations only for the L = 
300 h Mpc subvolumes (N cut = 8). The agreement with simulations is independent of 
subvolume size as long as the subvolume size is large enough for 5 to be in the linear regime, 
















46 


4. Measurement of position-dependent power spectrum 


and the window function is taken into account. 

We first compute the integrated bispectrum by using a model for the bispectrum in 


eq. (4.3) and perform the eight-dimensional integral. Because of the high dimensionality, 


we use the Monte Carlo integration routine in GNU Scientific Library to evaluate the 
angular-averaged integrated bispectrum. In the following, we consider two different models 
for the matter bispectrum. 


4.2.1 Standard perturbation theory 

The standard perturbation theory (SPT) [IT] gives the tree-level matter bispectrum as 

5 SPT (ki, k 2 , k 3 ) = 2[Pi(k 1 )Pi(k 2 )F 2 (k 1 , k 2 ) + 2 cyclic], (4.17) 


where Pi(k) is the linear matter power spectrum, and 


F 2 (ki k 2 ) = - + - kl ' k ' 2 
21 1? 2) 7 + 2 A h k 2 


ki h 
k 2 h 


2 

+ 7 


ki • k 2 

k\k 2 


(4.18) 


In order to normalize the integrated bispectrum, we need an expression for the mean 
subvolume power spectrum Pi,{k). For this we use the linear power spectrum convolved 
with the window function, 


P, 


i,M - Vl 


d 3 q 

(2tt) 3 


Pi(\k- q \)\W L (q)\^ 


(4.19) 


while the variance of the mean density fluctuation in the subvolumes is given by eq. (4.5). 
Both quantities are calculated through Monte Carlo integration. 

We compare the normalized integrated bispectrum measured from the simulations with 
the SPT prediction in figure [473] (red lines). The SPT prediction is independent of redshift. 
This is because the linear power spectra at various redshifts are only different by the 
wavenumber-independent linear growth factor, D 2 (z). Therefore, the linear growth factor 
cancels out in the normalized integrated bispectrum. The SPT predictions agree with 
the simulations relatively well at z > 1 and k < 0.2 h Mpc -1 , whereas they fail at lower 
redshifts as well as on smaller scales, where nonlinearities become too strong to be described 
by SPT. Especially, the BAO amplitudes at k > 0.2 h Mpc -1 are affected: while the SPT 
predictions are redshift-independent, the simulations show smaller BAO amplitudes at 
lower redshifts. 


The eight-dimensional integral in eq. (4.3) simplifies greatly if we focus on the squeezed- 


limit bispectrum. In appendix [Aj we show (note that I?spt = £>sqi 


d 2 Qi 


B S pt{ k - qi, -k + q 3 + q 3 , -q 3 ) 
Pl(k)Pi(q 3 ) + 0 


4:71 

68 1 din k 3 Pi(k) 

21 3 din A; 


gl,3V 
k ) 


(4.20) 
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Figure 4.3: The SPT and the F| ff (ki, k 2 ) predictions for the normalized integrated bispec¬ 
trum at different redshifts. The red solid and blue dot-dashed lines are computed by the 
direct integration of the eight-dimensional integral (eq. (4.3)) with the standard F 2 (ki, k 2 ) 
kernel and the linear power spectrum, and F% s (k l5 k 2 ) and the nonlinear power spectrum, 
respectively. The green dashed lines show the squeezed-limit approximation (eq. (4.21)) to 
the SPT results. The IV-body simulation results are shown by the black data points with 
the error bars showing the standard deviation on the mean measured from 160 simulations. 


for k q \, f/ 3 . We can then apply eq. (4.4) and perform all the integrals analytically in 
the limit of kL —> oo to obtain 


iB LtS p T (k) — —2 
L • 


d 2 a* 


d 3 q 1 f d 3 q 3 


47T J (27r) 3 J 


B S pt{ k - qi, -k + qi + q 3 , -q 3 ) 


x W L (c\i)W L (—c\i “q 3 )lpL(q 3 ) 


/cL—>• oo 


68 1 din k 3 Pi(k) 

21 3 din A; 


Pi(k)a 2 


L ■ 


(4.21) 


Comparing this result with eq. (4.7), we find that the linear response of the power spectrum 
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4. Measurement of position-dependent power spectrum 


in SPT is given by 


d In Pi(k) 


dd 


SPT 


68 

21 


1 din k 3 Pi(k) 
3 dink 


(4.22) 


The green dashed lines in figure |4.3 show the squeezed-limit approximation given in 
eq. (4.21). While they are different from the full integration (red solid lines) at k < 
0.2 h Mpc -1 , for which the squeezed-limit approximation fails and the direct integration 
is required, they agree well, with the fractional difference being less than 1.5% (1% for 
L = 600 h -1 Mpc), at k > 0.2 h Mpc -1 , corresponding to a value of l/{kL) < 0.02. Thus, 
the squeezed-limit is reached already with good precision for kL > 50. 


Eq. (4.21) does not contain any window function effect apart from that in the vari¬ 


ance a2 - While this is a good approximation for the slowly-varying part of the integrated 
bispectrum, it does not capture the smearing of the BAO features due to the window func¬ 
tion. We incorporate this effect by replacing din Pi{k)/dink with appropriately convolved 


forms, conv [dPi(k) / d In k] / conv[P/(fc)], in eq. (4.21). This form is motivated by the sepa¬ 


rate universe approach discussed in section 4.3 


and provides an accurate result as shown 


in figure 4.3 


4.2.2 Bispectrum fitting formula 

The SPT predictions fail on smaller scales as well as at lower redshifts where nonlinearity 
becomes too strong to be described by SPT. An empirical fitting formula for nonlinear 
evolution of the matter bispectrum was proposed in [ 140 j and further improved in [63]. In 
short, the form is the same as the tree-level matter bispectrum, but i^kijk^) is replaced 
by an effective kernel, iTf® (ki, k 2 ), which contains nine fitting parameters, {cq, • • • ag}, to 


account for nonlinearity (see eqs. 2.6 and 2.12 in [64] for details). Therefore, we use 
P| ff (k 1 ,k 2 ) and compute the integrated bispectrum by performing the eight-dimensional 
integral numerically with Monte Carlo integration. We use the same values of the best- 
fit parameters provided in table 2 of [64], which were calibrated by fitting to simulation 
results between z = 0 and z = 1.5. In contrast to the SPT formalism that uses the linear 


power spectrum in eq. (4.17), the fitting formula uses the nonlinear power spectrum, for 
which we use the mean power spectrum measured from the 160 simulation boxes. For 
the normalization of the integrated bispectrum, we convolve the nonlinear power spectrum 
with the subvolume window function as in eq. (4.19). Note that the fitting formula 


is not specifically designed for the squeezed configuration, but instead was calibrated to a 
wide range of triangle configurations of the matter bispectrum. 


The blue dot-dashed lines in figure |4.3| show the normalized integrated bispectrum 
computed with F 2 ff , which clearly depends on redshift. At z > 1, the modeling 
and the simulations are in good agreement at k < 0.2 h Mpc -1 . At k > 0.2 h Mpc -1 , 
although the F 2 eff modeling predicts larger broad-band power of the normalized integrated 
bispectrum, the BAO amplitudes still agree well with the simulations. This is most obvious 
for the two BAO peaks at 0.25 h Mpc -1 < k < 0.35 h Mpc -1 . On the other hand, at 
z = 0, the modeling predicts much larger normalized integrated bispectrum on small 
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scales than measured in the simulations, so that the fitting formula does not perform much 
better than tree-level perturbation theory ah z = 0. 


4.3 Separate universe approach 

In the second approach, we compute the effects of a long-wavelength density fluctuation on 
the small-scale power spectrum by treating each over- and under dense region as a separate 
universe with a different background density. This approach thus neglects the finite size of 
the subvolumes and is valid for wavenumbers which satisfy kL 1 (specifically, kL > 50 
for percent-level accuracy). 

The power spectrum in a separate universe with an infinite-wavelength density per¬ 
turbation, 5, with respect to the global flat ACDM cosmology can be expanded as in 
eq. (4.6). Through eqs. (4.7)—(4.8), the normalized integrated bispectrum is equal to the 


linear response of the nonlinear matter power spectrum at wavenumber k to 5: 

iB L (k) d\nP(k) 


P(k)aj 


d5 


(4.23) 


This is not exactly true if the subvolumes for which iB^k) is measured are not spheri¬ 
cal. For example, since the cubic window function is anisotropic, the integrated bispec¬ 
trum might pick up contributions from the tidal held. However, we have verified that the 
anisotropy of the cubic window function has a negligible effect, by computing the dipole 


and quadrupole of the integrated bispectrum through eq. (4.3). The ratios to the monopole 
are less than 1CT 5 on the scales of interest. 

A universe with an infinite-wavelength density perturbation with respect to a hat fidu¬ 
cial cosmology is equivalent to a universe with non-zero curvature. This alters the scale 
factor, Hubble rate, and linear growth as shown in chapter [3j and thus affects the power 
spectrum. Say this long-wavelength overdensity is 


m = 


p(t) 

pit) 


i = 


m 

Pit o) 


S(t 0 ) 


(4.24) 


where p(t) is the background matter density in the fiducial cosmology, D(t) is the linear 
growth factor in the same cosmology, p is the background matter density in a slightly curved 
universe, t 0 is a reference time such that a(t 0 ) = 1, and <5 0 is the density perturbation at 
t 0 . Note that as in chapter |3j we denote the quantities in the modified (curved) cosmology 
with a tilde. 


In eq. (4.24) and in the following of this section, we assume that 5(t) is small and evolves 
linearly. This is justified because here we consider the sub volume to be 300 h~ 1 
9x 10" 4 . 


a 


Mpc, and so 

L (z = 0) ~ 9 x 10 -4 . One can also see this in figure 4.1[ |h(r L )| < 0.08 in 512 subvolumes 
with L = 300 h~ l Mpc of one realization at z — 0. Therefore, we shall consider the effect 
on the power spectrum only to the linear order in S(t), and drop S n {t) for n > 2. 
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4. Measurement of position-dependent power spectrum 


With this assumption, the scale factor in the modified cosmology is given by 


d(t) 


a(t ) 



(4.25) 


Since the physical coordinates are the same in two cosmologies, this implies that the 
comoving coordinates of the two cosmologies are related by 


x = 


a(t) 

h(t) X = 


i + I'm 


x . 


(4.26) 


What we want is to compare the observables between the fiducial and modified cosmologies, 
so the quantities computed in the comoving coordinates of the modified cosmology have 
to be mapped to that with respect to the comoving coordinates of the fiducial cosmology. 
Specifically, in order to match to the comoving coordinates of the fiducial cosmology, we 
transform the comoving coordinates in the modified cosmology as 


x = 


1 - 3 m 


X = cx 


(4.27) 


where c is a constant at a given time. This assures x = x. 

Let us now consider how the transform in the coordinates affects the two-point statistics. 
Since the correlation function is a dimensionless scalar quantity, in the new coordinates 
x with the transformation of X{ = ryay (for i — 0, 1 , and 2 ) the correlation function £(x) 
must be 

f (#o, xi, x 2 ) = f (rr 0 , x u x 2 ) • (4.28) 

The power spectrum is then transformed as 

P(k 0 ,k u k 2 ) = j<Pxl(x 0 ,x 1 ,x 2 )e- i ^ +ilkl+i ^ ) 

= C 0 C 1 C 2 J d 3 x ^(a:o,a;i,a; 2 )e" i(coXoI ' 0+Cl3:i ^ 1+C2a;2 ^ ) 

= c oCl c 2 J d 3 x f( x 0 ,x 1: x 2 )e-^ X0ko+xlkl+X2k2) 

= c 0 CiC 2 P(k 0 , fci, k 2 ) = c 0 CiC 2 P(cofco, Cik^p) , (4.29) 


where we define p = cj l ki- This makes sense because k ~ and so k t ~ x^ 1 = c^ 1 x^ 1 ~ 
c7 t 1 kj (see also appendix A of [123]). 

Inserting c through eq. (4.27), we have the change in power spectrum up to the linear 
order of 5{t) as 


P(k,t) ->• 



r i_ /x i 

3 ( 

r i-i 

->■ 


Plk 
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I ICO 
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T -—1 


3 d In k 


= [l-S(i)] P(M) 

= t, ldlnep(k,t) 


1 din P(k, t ) 
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i— 
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din k 
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(4.30) 
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Eq. (4.30) is also known as the “dilation” effect in |97l 198]. which is the consequence that the 
presence of the long-wavelength overdensity perturbation slows down the local expansion 
rate. 

Another effect arises due to the change in the “reference density”. That is, the back¬ 
ground density in two cosmologies are related by pit) = p(t) [l + 5(t)], and since power 
spectrum is proportional to density squared, in the overdense universe the power spectrum 
at the linear order of 5{t) becomes 


P(k,t) -> [l + <S(()] 2 P(M) = [1 +2.5(f)] P(k,t) . 


(4.31) 


Combining the effects of dilation and the reference density, and using the scale factor 
instead of time, the power spectrum in the presence of 5 is given by 


P(k,a\5) = [1 + 2 6(t)]P(k,a) 


, ld]nk 3 P(k,t) jM 
1 ‘ 3 d\nk Ht) 


= P ( k, a 


1 +) 


1 + I 2 3 


1 din k 3 P(k, a) 


dink 


5(a) 


(4.32) 


Note that this expression is only valid to linear order in 5. 

Both P(k) and 5 are measured in a finite volume, described by the window function 
Wl- hi order to take this into account, eq. (4.32) is convolved by the window function. 


Note that we take the convolution after applying the derivative din k 3 P(k)/dink, rather 
than taking the derivative of the convolved power spectrum. This is because the window 
function is fixed in terms of observed coordinates (in the fiducial cosmology), i.e., it is not 


subject to the rescaling of eq. (4.26). Taking the slope of the convolved power spectrum 
would correspond to a window function defined in the “local” curved cosmology. 


4.3.1 Linear power spectrum 

For the linear power spectrum, Pi, we have 


Pi 



D(a\l- §£(«)]) 
D(a) 


Pi{k,a ) 


(4.33) 


As shown in section 3.2.2 (see also appendix D of HU), the linear growth factor is changed 
following 


D I a 


1 - -6{a) 


= D(a ) 


1 + —5(a) 


(4.34) 


where D{a) is the growth factor in the fiducial cosmology. The prefactor 13/21 is only 
strictly valid for an Einstcin-de Sitter cosmology; however, the cosmology dependence 
is very mild. The fractional difference of dlnD{a)/d8 between ACDM cosmology and 
Einstein-de Sitter universe at z = 0 is at the 0.1% level. 
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4. Measurement of position-dependent power spectrum 


Putting everything together, eq. (4.32) yields for the linear response function of the 
linear power spectrum 


d In Pi(k, a) 68 1 din k 3 Pi(k, a) 

” 21 _ 3 


(4.35) 


dS(a) 21 3 dink 

This result (which again is only exact for Einstein-de Sitter) matches the expression derived 


from the F 2 kernel given in eq. (4.22). 


4.3.2 SPT 1-loop power spectrum 

Expanding matter density fluctuations to third order, one obtains the so-called “SPT 1- 
loop power spectrum” given by P(k, a) = Pi(k, a) + P 22 (k, a) + 2 P ] 3 (/c, a), where [Tf] 


P-n(k,a) = 2 / fl(9. o)fl(|k - q|, o) [F 2 (q, k - q)] 2 , 


( 2*0 

27T k 2 f c 

2P W (k,a) = —P l (k,a) J 


(4.36) 


dq 

W> 


P L (q, a) 


X 


m l _ 158+ 12* - 42 d + m - V) 1. (j|^ 


Both P 22 and P 13 are proportional to D 4 (a). Modifying the growth factor as described in 
section 4.3.1[ we obtain the linear response function of the SPT 1-loop power spectrum as 


d\nP(k,a ) 68 


dS(a) 


21 


1 din k 3 P(k, a ) 26 P 22 {k , a) + 2P\ :i (k, a) 

3 din A: ”*”21 P(k,a ) 


(4.37) 


Note that this can easily be generalized to n loops in perturbation theory by using that 
dlnP( ra _i oop )(A;, a)/d\nD(a) = 2n + 2. We include the window function effect by computing 
conv [dP(k) / dS\ / conv[P( k)\. 

Figure |4.4| compares the linear theory and the SPT 1-loop predictions with the N- 
body simulation results. The SPT 1-loop prediction captures the damping of BAOs due 
to nonlinear evolution, and agrees well with the simulation results at z — 1, 2, and 3. 
This is expected from the excellent performance of the 1-loop matter power spectrum at 
high redshifts as demonstrated by m The agreement degrades rapidly at z — 0, also as 
expected. Note that comparing z — 2 and 3, the 1-loop prediction seems to agree better 


with the measurements at z = 2. However, as mentioned in section 4.1, transients and 


other systematics might have an impact of up to 5% on the measurements at z = 3, which 


is larger than the difference shown in the top left panel of figure 4.4 


4.3.3 halofit and Coyote emulator 

We now apply the separate universe approach to simulation-calibrated fitting formulae for 
the nonlinear matter power spectrum, specifically the halofit prescription [ 152 ] and the 
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Figure 4.4: Normalized integrated bispectrum from the IV-body simulations (points with 
error bars) and the linear response functions, dlnP(k, a)/d5(a), computed from the sep¬ 
arate universe approach combined with perturbation theory. The red dashed lines show 


the linear theory results (eq. (4.35)), while the blue solid lines show the SPT 1-loop results 
(eq. (4.37)). The agreement between the 1-loop predictions and the simulation results is 


very good at z > 1. Note that the difference between the normalized integrated bispectrum 
and the linear response function at k < 0.2 h Mpc^ 1 is due to the squeezed limit not being 
reached yet (see the text below eq. (4.22[)). 


Coyote emulator im These prescriptions yield P(fc, a) for a given set of cosmological 


parameters, so that eq. (4.32) can be immediately applied. However, the Coyote emulator 


does not provide predictions for curved cosmologies, and we hence adopt a simpler approach 
here. 

In case of the linear power spectrum, the effect of the modified cosmology enters only 


through the modified growth factor given in eq. (4.34). Correspondingly, we can approx¬ 


imate the effect on the nonlinear power spectrum by a change in the value of the power 
spectrum normalization erg at redshift zero, 


' 13 j 

1+ 2l' 50 


08 , 


cr 8 —y 


(4.38) 
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4. Measurement of position-dependent power spectrum 



Figure 4.5: Same as figure [F4| but for the linear response functions computed from halofit 
(red solid), the Coyote emulator (green dat-dashed), and the halo model (blue dashed). 


where we have used the Einstein-de Sitter prediction. Therefore, the nonlinear power 
spectrum response becomes 


d In P n i ( k , a) 13 din P n i(k, a) 


d5(a) 


21 d In erg 


2 _ 1 d\nk 3 P n i(k, a) 
3 d In k 


(4.39) 


The results of applying eq. (4.39) to halofit (red solid) and the Coyote emulator (green 
dot-dashed) are shown in figure 4.5 In terms of broad-band power, the halofit prediction 
provides a good match. However, the predicted BAO amplitude are larger than the mea¬ 
surement, especially at low redshift at k > 0.3 h Mpc -1 . Also, while the BAO phases 
of halofit follow the SPT prediction, there are some differences with respect to the mea¬ 
surement of the IV-body simulations due to the nonlinear evolution. The Coyote emulator 
performs to better than ~ 2% over the entire range of scales and redshifts. It slightly 
underpredicts the small-scale power at k > 0.3 h Mpc -1 for z > 1. For redshifts z > 2 and 
on the scales considered, the 1-loop predictions are of comparable accuracy to the Coyote 
emulator, while the latter provides a better fit at lower redshifts. Finally, note also our 


previous caveat regarding transients at the end of section 4.3.2 
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4.3.4 Halo model 

In the halo model (see [35] for a review), all matter is assumed to be contained within 
halos with a certain distribution of mass given by the mass function, and a certain density 
profile. Along with the clustering properties of the halos, these quantities then determine 
the statistics of the matter density held on all scales including the nonlinear regime. Ap¬ 
point functions can be conveniently decomposed into one- through TV-halo pieces. In the 
following, we will follow the most common halo model approach and assume a linear local 
bias of the halos. 

Adopting the notation of [157] . the halo model power spectrum, P HM (/c), is given by 
P HM (k) = P 2h (k) + P lh (k) , P 2h (k) = [Il(k)] 2 P t {k) , P lh (k) = I$(k,k) , (4.40) 

where 

/”(fci,-• • ,k m ) = /dlnMn(lnM)(|) b n (M)u(M\ki) ■ ■ ■ u(M\k m ) , (4.41) 


and rt(ln M) is the mass function (comoving number density per interval in log mass), 
M is the halo mass, b n (M) is the n-th order local bias parameter, and u(M\k) is the 
dimensionless Fourier transform of the halo density profile, for which we use the NFW 
profile [ 118] . We normalize u so that u(M\k —» 0) = 1. The notation given in eq. (4.41) 
assumes b 0 = 1. u(M\k) depends on M through the scale radius r s , which in turn is 
given through the mass-concentration relation. All functions of M in eq. (4.41) are also 
functions of z although we have not shown this for clarity. In the following, we adopt the 
Sheth-Tormen mass function [135] with the corresponding peak-background split bias, and 
the mass-concentration relation of |24j . The exact choice of the latter has negligible impact 
on the mildly nonlinear scales, but does not affect the conclusion. 

We now derive how the power spectrum given in eq. (4.40) responds to an infinitely 
long-wavelength density perturbation 5, as was done for the haloht and Coyote emulator 
approaches. For this, we consider the one-halo and two-halo terms separately. The key 
physical assumption we make is that halo profiles in physical coordinates are unchanged 
by the long-wavelength density perturbation. That is, halos at a given mass M in the 
presence of 5 have the same scale radius r s and scale density p(r s ) as in the fiducial 
cosmology. This assumption, which is related to the stable clustering hypothesis, can be 
tested independently with simulations, but we shall leave it for future work. Given this 
assumption, the density perturbation 5 then mainly affects the linear power spectrum, 
which determines the halo-halo clustering (two-halo term), and the abundance of halos at 
a given mass. 

We begin with the two-halo term. The response of the linear power spectrum is given 
by eq. (4.35). The expression for the two-halo term in eq. (4.40) is simply the convolution 
(in configuration space) of the halo correlation function in the linear bias model with 
the halo density profiles. By assumption, the density profiles do not change, hence l\ 
only changes through the bias b\{M ) and the mass function n(lnM). The bias bjy(M) 
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quantifies the iV-th order response of the mass function n(lnM) to an infinite-wavelength 
density perturbation mm-- 


b N (M) = 


d N n( In M) 


ra(lnAf) dS N 


<5=0 


(4.42) 


We then have 


dn{\n M ) 


85 

<%i (Af) 


c>5 


<5=0 


<5=0 


= b 1 (M)nQnM) , 

-[6!(M)] 2 + 6 2 (M) . 


(4.43) 


Thus, 


^ll(k) = I dlnMn(hiM) |[ &1 (M)] 2 - [^(M)] 2 + 6 2 (M)} u(M|fc) 

= j dlnAf ra(lnM) 6 2 (Af)u(M|fc) 

= / 2 (fc). (4.44) 

In the large-scale limit, k —y 0, this vanishes by way of the halo model consistency relation 


d\nMn(\nM)(jjb N {M) = { J’ ^ > j ’ 


(4.45) 


For finite k however, eq. (4.44) does not vanish. Thus, the linear response function of the 
two-halo term becomes 


P 2l \k) + 2ll(k)l\(k)Pi(k) . (4.46) 


dp 2h (k) 


68 1 din k 3 Pi(k) 

d5 

<5=0 

21 3 dink 


Note that we recover the tree-level result given in eq. (4.35) in the large-scale limit. 
Strictly speaking, this expression is not consistent, since the term I 2 implies a non-zero 


bo while in eq. (4.40) we have assumed a pure linear bias. Of course, if we allowed for 
b 2 in eq. (4.40), we would obtain a contribution from b 3 in eq. (4.46), and so on. This 
reflects the fact that the halo model itself cannot be made entirely self-consistent. Note 


that in eq. (4.46) the slope is taken from the linear, not two-halo power spectrum. This is 


a consequence of our assumption that halo profiles do not change due to 5 ; in other words, 
having din k 3 P 2h /d\nk would imply that the profiles do change (in the sense that they are 
fixed in comoving, rather than physical coordinates). 

We now turn to the one-halo term. Given our assumption about density profiles, this 
term is much simpler. The only effect is the change in the mass function, which through 


eq. (4.42) (for N — 1) yields 


~I° 2 (k,k) = Ii(k,k) . 


( 4 . 47 ) 
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We thus obtain 


dP lh (k) 


dd 


Putting everything together, we obtain 


< 5=0 


= ll(k,k). 


(4.48) 


dlnP HM (k) 


dd 


( 5=0 


= [P RM (k)] 


-i 


68 

21 


1 din 


k 3 Pi(k) \ 


dink ) 


The prediction of eq. (4.49) is shown as the blue dashed lines in figure 4.5 


P 2h (k) + 2lf(k)I 1 1 (k)P l (k ) + l\{k, k) 

(4.49) 
The amplitude 


and broad-band shape agree with the simulations well. The main discrepancy in the halo 
model prediction is the insufficient damping of the BAO wiggles. 

An alternative approach to derive the halo model prediction for iBi(k ) is to use higher 
iV-point functions [Ml 97] . which are decomposed into one-, ..., N —halo terms. We now 
compare eq. (4.49) with the results of pT7J, which were derived from the halo model four- 


point function in the collapsed limit, 
approaches. Their eq. (27) is 


Note that the squeezed limit is assumed in both 


dlnP HM (k) 


dd 


5=0 


= [ P RM (k )] 


-1 


21 


idink 3 p 2 py>\ k ) + ii(k, k ) 


dink 


(4.50) 


There are two differences to eq. (4.49): the term oc If is absent, and the slope is taken 
from from P 2h rather than P[. The If term is absent in eq. (4.50) as by assumption b 2 
was taken to be zero in the four-point function of |9TJ; as discussed above, its inclusion is 
somewhat ambiguous given the lack of self-consistency of the halo model approach. The 
different power spectrum slopes are due to the different sources of this term in the two 
derivations. In our case, the assumption of unchanged halo profiles dictates the form of 


eq. (4.49). In the derivation of eq. (4.50), the slope originates from the integral over the F 2 
kernel in the three-halo term, which proceeds as described in appendi x [A| b ut involves P 2h 
instead of Pi. Note however that the numerical difference between eq. (4.50) and eq. (4.49) 
is only at the percent level. 


4.4 Dependence on cosmological parameters 

Both the matter power spectrum and (integrated) bispectrum depend on the cosmological 
parameters such as Q m , cr 8 , n s - However, the normalized integrated bispectrum is much less 
sensitive to cosmology as the leading cosmology dependence is taken out by the normalizing 
denominator. 


Eq. (4.37) is useful for understanding the dependence of the response function of the 


power spectrum (and thus the normalized integrated bispectrum) on cosmological param¬ 
eters. The second term depends on the local spectral index of the matter power spectrum, 
din k 3 P[k)/dink, which depends on the initial power spectrum tilt, n s , and the matter 
and radiation densities which change the redshift of matter-radiation equality as well as 
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4. Measurement of position-dependent power spectrum 



Figure 4.6: The linear response functions computed from the SPT 1-loop power spectrum 
with various cosmological parameters at z — 2. The fiducial cosmology (fl m = 0.27, 
as = 0.7913, and n s = 0.95) is shown in green solid lines. The red dot-dashed (blue 
dashed) lines represent the cosmologies with -5% (+5%) of the fiducial parameters, Q m 
(left), erg (middle), and n s (right). 


the BAO scale. It also depends on the shape of BAO wiggles, and increasing the ampli¬ 
tude of the matter power spectrum (<r 8 ) leads to a stronger damping of the BAO feature. 
Increasing <jg further increases the last term, which is proportional to cr|. 

Figure 4.6 shows the linear response functions, dha.P(k,a)/dS(a) computed from the 
SPT 1-loop power spectrum (eq. (4.37)) at z = 2 when varying cosmological parameters 
by ±5%. The effects on the response functions are at the percent level or less, illustrating 
the weak cosmology dependence of this observable. On the scales considered, the shift in 
the BAO scale when varying Q m leads to the relatively largest effect. We expect that the 
sensitivity to changes in <jg will be higher on smaller, more nonlinear scales. 


4.5 Fisher matrix calculation 


Now that we understand the behavior of the integrated bispectrum, how does it compare 
with the full bispectrum in terms of measuring the cosmological parameters, particularly 
the primordial non-Gaussianity? In this section, we perform the Fisher matrix calculation 
(see e.g. [ISP] for a review) for both the full bispectrum analysis and the integrated 
bispectrum technique, and compare the performances between the two methods. We shall 
use the simplest primordial non-Gaussianity model as discussed in section 2.1.4, i.e. 


B g { ki, k 2 , k 3 ) — &i-BspT(ki, k 2 , k 3 ) + blb 2 Bb 2 (ki, k 2 , k 3 ) + &i/NL-B/ NL (ki, k 2 , k 3 ) 

iB L}9 (k) = bliB L £ PT (k) + b\b 2 iB LM {k) + blf NL iB L j NIj (k) . (4.51) 


The Fisher matrix of the reduced bispectrum Q(k\, k 2 , k 3 ) = B(k\, k 2 , k 3 )/[P{ki)P[k 2 )+ 
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2 cyclic] is given by 


F Q,<*P ~ ^ 

ki,k2,k 3 <k 

n 


dQ(k i, k 2 , fa) dQ(fa, fa, fa) 


dp 0 


dpp AQ 2 (fa, fa, fa) 


(4.52) 


where (fa, fa, fa) have to form a triangle, p a G [5i,&2,/nl] are the parameters we want 
to constrain, and AQ 2 (fa,fa,fa) is the variance of the reduced bispectrum estimator. 
Similarly, the Fisher matrix of the integrated bispectrum can be written as 


Fi, 


ibL,ct(3 


E E 


dibfak) dib L (k) 1 
dp a dpf Aib 2 L (k) 


(4.53) 


Here, we ignore off-diagonal elements of the covariance matrix of Q or ibi- In general, 
nonlinear evolution generates non-vanishing covariances. This is, however, justified at high 


redshift (z > 2) as we show in the right panel of figure B.l 


The variance of the estimator for the integrated bispectrum is computed in appendix [Bj 
The variance of the estimator for the bispectrum is computed [144j 


AB 2 (fa, fa, fa) = 


7TS123 

fa fa fa 


P(fa)P(fa)P(fa) 


(4.54) 


where S 123 = 6 , 2 , 1 for equilateral, isosceles, and general triangles, respectively (we set A k 
to be the fundamental frequency). Similar to the integrated bispectrum, we assume that 
the variance of the reduced bispectrum is dominated by the numerators, so 


A Q 2 (fa,fa,fa) 


AB 2 (fa,fa,fa) 
[P(fa)P(fa) + 2 cyclic ] 2 


(4.55) 


Since galaxies are observed in redshift space, we model the redshift-space distortions by 
the simple Kaiser factor, P z = K p P r and B z = K b B r , where the subscripts r and z denote 
the real- and redshift-space quantities, and 


2 / 1 
An = 1 + + - 

p 3 fa 5 



K b 


2 / 

1 + -— + 
3 61 


1 

9 



(4.56) 


with / = din D/ din a being the growth rate. The derivatives of Q and ib^ with respect to 
b\ thus contain the contributions from dK p /db\ and dK b /db\. The variances of the redshift- 
space reduced bispectrum and the normalized integrated bispectrum with the Poisson shot 
noise are then given by 


AQ 2 (fa, fa, fa) 


Aib 2 L (k) 


KS 123 [Pz(k l) + Pshot][Pz(fa) + Pshot][P Z (fa) + -P s hot] 
fafafa [P z (fa)P z (fa) + 2 cyclic ] 2 

Vl K, z + P sho t/V L ][P L , z (k) + -Ashot] 


a 


\, z Plz(k) 


V r N kL 


(4.57) 
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Figure 4.7: Two-dimensional joint 95% C.L. constraints on galaxy bias and primordial non- 
Gaussianity for BOSS (left) and HETDEX (right). The survey parameters are in the top 
right panel. The top left, bottom left, and bottom right panels show the joint constraints 
on (6i, b 2 ), (b\, /nl), and (b 2 , /nl) marginalized over /nl, b 2 , and bi, respectively. The blue 
dashed, green solid, and red dot-dashed lines are for full bispectrum with /c min — kp — 
27t/ 14- 1//3 , full bispectrum with k min = kp,L = 27 t/L where L is the largest subvolume size 
(600 h _1 Mpc), and integrated bispectrum for six sizes of subvolumes (100 h -1 Mpc to 
600 h -1 Mpc with an increment of 100 h~ l Mpc), respectively. 

where P z (k) = b\K p Pi(k ), a\ z = b\K p a 2 L , and P L , z (k) = b\K p P L (k). 

Figure E3 shows the two-dimensional joint 95% C.L. constraints on galaxy bias and 
primordial non-Gaussianity for BOSS [4] (left) and HETDEX [71] (right). The survey 
parameters and the fiducial cosmological parameters are shown in the top-right of each 
panel. The blue dashed line is the full bispectrum analysis with fc min = 2tt jV r 7 being the 
fundamental frequency of the entire survey V r \ the green solid line is the full bispectrum 
analysis with k m j n = 2it /V L being the fundamental frequency of the largest subvolume for 
the integrated bispectrum (Vp = [600 h~ l Mpc] 3 ); the red dot-dashed line is the integrated 
bispectrum with six sizes of subvolumes from 100 h~ l Mpc to 600 h~ l Mpc (with an 
increment of 100 h Mpc) with k m \ n = 2n/V^' i being the fundamental frequency of the 
corresponding subvolumes. 

One Ends that as long as k m[n is set to be the fundamental frequency of the largest 
subvolume, the integrated bispectrum technique gives the similar constraint on /nl com¬ 
pared to the full bispectrum analysis^ On the other hand, the integrated bispectrum has 

2 Note that the numerical results are sensitive to the choices of k m \ n and fc max because we count the 

Fourier inodes in this range. For different lines, although fc max is set to be the same, in practice we stop 
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poor constraints on b\ and b 2 compared to the full bispectrum analysis. In particular, in 
the top-left panel the signal of integrated bispectrum has a strong degeneracy between 
b\ and b 2 . This is somewhat expected because in figure |2.1 one finds that not only the 
black solid and black dashed lines have similar scale dependences in a given subvolume, 
but in different subvolumes they all have similar contribution, unlike the bispectrum of the 
primordial non-Gaussianity. This thus makes it difficult to break the degeneracy between 
61 and b 2 using the integrated bispectrum. 

We also find that while the full bispectrum analysis and the integrated bispectrum 
technique give similar constraints on /nl, the number of counted Fourier modes differ 
dramatically. For example, for the BOSS parameter, the full bispectrum analysis with 
kmin = 27t/ 14- 1 ^ 3 counts 7113 configurations of triangles, whereas the integrated bispectrum 
technique counts only 54 Fourier modes. Even if k m ; n is set to be 2tt/V l ' for the full 
bispectrum analysis, there are still 6730 configurations of triangles. This is a big advantage 
for the integrated bispectrum technique because estimating the covariance matrix from 
mock catalogs would require a large number of realizations, which will be difficult (but not 
impossible) to obtain for the full bispectrum. This difference in the number of counted 
Fourier modes also explains why the full bispectrum analysis has much better constraints 
on bi and b 2 . However, many of the triangles do not contain much more information on 
/nl, so if one is interested in measuring / NL , the integrated bispectrum technique provides 
an easier approach and captures most of the information. 


4.6 Discussion and conclusion 


In this chapter, we have demonstrated a novel method to measure the squeezed-limit bis¬ 
pectrum. By measuring the correlation between the mean density fluctuation and the 
position-dependent power spectrum, we obtain a measurement of an integral of the bispec¬ 
trum (integrated bispectrum) without having to actually measure three-point correlations 
in the data. The integrated bispectrum is dominated by the squeezed-limit bispectrum, 
which is much easier to model than the full bispectrum for all configurations. This is 
evidenced by figure |4~4| and figure 473, where we show model predictions accurate to a few 
percent using existing techniques and without tuning any parameters. 

A further, key advantage of this new observable is that both the mean density fluctu¬ 
ation and the power spectrum are significantly easier to measure in actual surveys than 
the bispectrum in terms of survey selection functions. I 11 particular, the procedures de¬ 
veloped for power spectrum estimation can be directly applied to the measurement of the 
position-dependent power spectrum. Additionally, the position-dependent power spectrum 
depends on only one wavenumber (at fixed size of the subvolume) rather than the three 
wavenumbers of the bispectrum. Consequently, the covariance matrix also becomes easier 


counting Fourier modes if k > fc max . Therefore they have different “true” fc max , and the contour area 
would be affected, especially the green solid lines seem to have slightly worse constraints on /nl compared 
to the red dot-dashed lines. Here, however, what we are interested in is the general properties between 
different methods, so we should neglect the effect due to different fc max . 
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4. Measurement of position-dependent power spectrum 


to model. 

We have measured the position-dependent power spectrum in 160 collisionless IV-body 
simulations with Gaussian initial conditions, and have used two different approaches, bis¬ 
pectrum modeling and the separate universe approach, to model the measurements. All of 
the approaches work well on large scales, k < 0.2 h Mpc^ 1 , and at high redshift. On small 
scales, where nonlinearities become important, the separate universe approach (section [4.3[ ) 
applied through the Coyote emulator prescription performs best at redshifts z < 2, while 
the SPT 1-loop predictions perform equally well at z > 2. Both show agreement to within 
a few percent up to k — 0.4 h Mpc -1 . Accurate predictions for the position-dependent 
power spectrum on these and even smaller scales can be obtained by applying the separate 
universe approach to dedicated small-box IV-body simulations of curved cosmologies, as 
described in section 16.21 

The normalized integrated bispectrum is relatively insensitive to changes in cosmo¬ 
logical parameters (section 4.4), and we do not expect that it will allow for competitive 
cosmology constraints. On the other hand, this property can also be an advantage: since 
this observable can be predicted accurately without requiring a precise knowledge of the 
cosmology, it can serve as a useful systematics test for example in weak lensing surveys. As 
an example, consider eq. (4.4) applied to shear measurements. A constant multiplicative 
bias 1 + m in the shear estimation contributes a factor (1 + m) 3 on the left hand side 
of the equation, and a factor (1 + m) 4 on the right hand side. Thus, by comparing the 
measured normalized integrated bispectrum with the (essentially cosmology-independent) 
expectation, one can constrain the multiplicative shear bias. 

The position-dependent power spectrum can also naturally be applied to the case of 
spectroscopic galaxy surveys, in which case the nonlinear bias of the observed tracers 
also contributes to the bispectrum and position-dependent power spectrum. Thus, when 
applied to halos or galaxies, this observable can serve as an independent probe of the 
bias parameters and break degeneracies between bias and growth which are present when 
only considering the halo or galaxy power spectrum. We shall exploit this to measure the 
nonlinear bias of the BOSS CMASS galaxies in chapter [5j 

We finally use the Fisher matrix to show that if multiple sizes of subvolumes are used, 
the position-dependent power spectrum captures most of the information of the local- 
type non-Gaussianity contained in the full bispectrum analysis, but for much less Fourier 
modes. This is a huge advantage for this novel technique as the computational requirement 
for estimating the covariance matrix is largely alleviated. Consequently, the position- 
dependent power spectrum provides an easier approach for hunting the primordial non- 
Gaussianity for future galaxy surveys. 








Chapter 5 

Measurement of position-dependent 
correlation function 


In this chapter, we measure the position-dependent correlation function and the integrated 
three-point function from real data, the SDSS-III Baryon Oscillation Spectroscopic Survey 
Data Release 10 (BOSS DR10) CMASS sample [3JE]. 

As introduced in section 2.2 the correlation between the position-dependent correlation 
function, 


l(r ’ ri) = / = S 

and the mean overdensity, 


d 3 x 5(r+x)h(x)IR L (r+x-r i )IR i (x-r i ) , (5.1) 


S(ri) = u./ 

is the integrated three-point function 


d 3 q 
(2vr) 3 


5(-q)W L (q)e 


-ir L q 


(5.2) 


Kl{t) = (£{r, t l )S(t l )) 

= ^2 f J J C(r + Xi + r L ,Xi + r L ,x 2 + r L ) 

x W L (r + x 1 )Wi(x 1 )Wi(x 2 ) 


(5.3) 


where Vl is the size of the subvolume. Inspired by the behavior in the squeezed-limit 
where r ^ L, we define the normalized integrated three-point function as Ar( r ') I G \ with 
a\ being the variance of the fluctuations in Vl- Note that when comparing the model to 
the measurements, we shall divide the model by /l, bndr y ( r ), which is given in eq. (2.37), to 
correct for the boundary effect. 

The integrated three-point function is simply the Fourier transform of the integrated 
bispectrum 

iB L (k)sinc(kr) . 


Kl(t) = 


2tt 2 


(5.4) 
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5. Measurement of position-dependent correlation function 


Therefore, if we do not have the analytical expression for the three-point function, we can 
first compute the integrated bispectrum and Fourier transform it to obtain the integrated 
three-point function. One example is the redshift-space integrated three-point function, for 
which we first evaluate the redshift-space integrated bispectrum with the explicit expression 
of the SPT redshift-space bispectrum given in appendix A.2 and then apply eq. (5.4). We 
also show that the precision of this operation (nine-dimensional integral in total) is within 
2% on the scales of interest (30 h Mpc < r < 78 /r _1 Mpc, which we will justify in 
section 5.1.3). 


This chapter is organized as follows. In section 5.1 and section 5.2, we measure the 


position-dependent correlation function from PTHalos mock catalogs [ 143 . 11041 HUB] and 
the BOSS DR10 CMASS sample, respectively. The cosmological interpretation of the 
measurement is in section 15.31 We conclude in section 15.41 


5.1 Measurement of PTHalos mock catalogs 


We first apply the position-dependent correlation function technique to the 600 PTHa¬ 
los mock galaxy catalogs of the BOSS DR10 CMASS sample in the North Galactic Cap 
(NGC). From now on, we refer to the real and mock BOSS DR10 CMASS samples as the 
“observations” and “mocks”, respectively. 

We use the redshift range of 0.43 < z < 0.7, and each realization of mocks contains 
roughly 400,000 galaxies. We convert the positions of galaxies in RA, DEC, and redshift 
to comoving distances using the cosmological parameters of the mocks. The mocks have 
the same observational conditions as the observations, and we correct the observational 
systematics by weighting each galaxy differently. Specifically, we upweight a galaxy if 
its nearest neighbor has a redshift failure (w z f) or a missing redshift due to a close pair 
(w cp ). We further apply weights to correct for the correlation between the number density 
of the observed galaxies and stellar density (w star ) and seeing (w see ). We apply the same 
weights as done in the analyses of the BOSS collaboration, namely FKP weighting, wfkp = 
[1 + P^h^jcomp] -1 [56], where P w = 20000 h~ 3 Mpc 3 , and n(z) and “comp” are the 
expected galaxy number density and the survey completeness, respectively, provided in 
the catalogs. Therefore, each galaxy is weighted by wross = ( w cp + w z f ~ l)'Wktar'Wkee w FKP- 


In this section, we present measurements from mocks in real space in section 5.1.3 and 


redshift space in section [5.1.4| The application to the CMASS DR10 sample is the subject 
of section 15.21 


5.1.1 Dividing the subvolumes 

We use SDSSPixQ to pixelize the DR10 survey area. In short, at the lowest resolution 
(res=l) SDSSPix divides the sphere equally into n x = 36 longitudinal slices across the 
hemisphere (at equator each slice is 10 degrees wide), and each slice is divided into n y = 13 
pieces along constant latitudes with equal area. Thus, for res=l there are n x x n y = 468 

1 SDSSPix: http://dls.physics.ucdavis.edu/~scranton/SDSSPix 
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RA [degree] RA [degree] 


Figure 5.1: Division of random samples into subvolumes with two resolutions in the RA- 
DEC plane. Each colored pattern extends over the redshift direction. 


pixels. In general the total number of pixels is n' x x n' y = (res n x ) x (res n y ) = (res) 2 x 468, 
and in this chapter we shall set res=1024. After the pixelization, the i th object (a galaxy 
or a random sample) has the pixel number (■ i x ,i y ). 

We use two different subvolume sizes. To cut the irregular survey volume into subvol¬ 
umes with roughly the same size, we first divide the random samples at all redshifts into 10 
and 20 slices across longitudes with similar numbers of random samples; we then divide the 
random samples in each slice into 5 and 10 segments across latitudes with similar numbers 
of random samples. Figure |5.1| shows the two resolutions of our subvolumes before the red- 
shift cuts. (Note that this resolution is different from the resolution of SDSSPix, which we 
always set to res=1024.) Each colored pattern extends over the redshift direction. Finally, 
we divide the two resolutions into three (z cut = 0.5108, 0.5717) and five (z cut = 0.48710, 
0.52235, 0.55825, 0.60435) redshift bins. 

As a result, there are 150 and 1000 subvolumes for the low and high resolution configu¬ 
rations, respectively. The sizes of the subvolumes are approximately V\j 3 = 220 h Mpc 
and 120 /r _1 Mpc, respectively^] The fractional differences between the numbers of the 
random samples in subvolumes for the low and high resolutions are within 1°'and 
-i 83 %’ respectively. Since the number of random samples represents the effective volume, 
all subvolumes at a given resolution have similar effective volumes. We assign galaxies into 
subvolumes following the division of random samples. 


2 The shapes of the subvolumes are not exactly cubes. For example, for the high resolution, the ratios 

of square root of the area to the depth, y/L x L y /L z , are roughly 0.78, 1.42, 1.51, 1.28, and 0.71, from the 
lowest to the highest redshift bins. The results are not sensitive to the exact shape of the subvolumes, as 
long as the separation of the position-dependent correlation function that we are interested in is sufficiently 
smaller than L x , L y , and L z . 
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5. Measurement of position-dependent correlation function 


5.1.2 Estimators in the subvolumes 


In the i th subvolume, we measure the mean overdensity with respect to the entire NGC, Si, 
and the position-dependent correlation function, £j(r). The mean overdensity is estimated 
by comparing the total weighted galaxies to the expected number density given by the 
random samples, i.e., 


Si 


1 Wg.i 

& ^r,i 


-i, 


viV s 

= E»= 1 Wg,j 
Ei= 1 W r,i 


W g,tot 
Wr, tot 


(5.5) 


where w g i and w Tji are the total weights (wboss) °f galaxies and random samples in the 
i th subvolume, respectively, and N s is the number of subvolumes. 

We use the Landy-Szalay estimator [91] to estimate the position-dependent correlation 
function as 


£ls ,i(r,n) = 


DDiir, n) ( Er^r,i] 2 -E 
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RRi(r,H ) T,n W g,i'Lr W r,i 


+ 1 


(5.6) 

where DDi(r,n), DRi(r,n), and RR t (r, //) are the weighted numbers of galaxy-galaxy, 
galaxy-random, and random-random pairs within the i th subvolume, respectively, and /i is 
the cosine between the line-of-sight vector and the vector connecting galaxy pairs (ri — r 2 ). 
The summations such as Wr > i an d E g w g,i denote the sum over all the random sam¬ 
ples and galaxies within the i th subvolume, respectively. The angular average correlation 
function is then £ LS ,i( r ) = f 0 (l l l £ls ,iir,n). 

Eq. (5.6) estimates the correlation function assuming that the density fluctuation is 
measured relative to the local mean. However, the position-dependent correlation function 
defined in eq. (5.1) uses the density fluctuation relative to the global mean. These two 


fluctuations can be related by <5 g i 0 bai = (1 + <5)<5iocai + S with <5 = ni oca i/n g i 0 b a i — 1- Thus, the 
position-dependent correlation function, £j(r), is related to the Landy-Szalay estimator as 


= (i + + S; 


(5.7) 


To compute the average quantities over all subvolumes, we weight by w r , t in the corre¬ 
sponding subvolume. For example, for a given variable g t in the i th subvolume, the average 
over all subvolumes, g , is defined by 
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^V, tot 


N s 

i=1 


(5.8) 


Since the number of random samples in each subvolume represents the effective volume, 
the average quantities are effective-volume weighted. Eq. (5.8) assures that the mean of 
the individual subvolume overdensities is zero, 
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We also confirm that £(r) from eq. (5.7) agrees with the two-point function of all galaxies 
in the entire survey, on scales smaller than the subvolume size. With £j(r) and <5,, we 
estimate the shot-noise-corrected integrated three-point function in the subvolume of size 
L as 


i((r) = 


1 


N a 

^V,tot , 

’ i=i 


ii(r)5i ~ 2f(r) 


(1 + a) E r W r,i 

Oi Yjt h r ,iCOmp rjUJ^i 
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n r iComp r 


w r 


(5.i°) 

where the second term in the parentheses is the shot noise contribution, and n Tt i and 
comp r i are the expected galaxy number density and the survey completeness, respectively, 
of the random samples. Similarly, we estimate the shot-noise-corrected variance of the 
fluctuations in the subvolumes of size L as 
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(5.11) 


where the second term in the parentheses is the shot noise contribution. We find that the 
shot noise is subdominant (less than 10%) in both i((r) and cr\. 


5.1.3 Measurements in real space 


Figure 5.2| shows the measurements of the two-point function £(r) from the entire survey 
(top left) and the normalized integrated three-point functions (bottom panels), 
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(5.12) 


for the subvolumes of two sizes (220 h Mpc in the bottom-left and 120 h _1 Mpc in the 
bottom-right panels). The gray lines show individual realizations, while the dashed lines 
show the mean. 

We now fit models of £(r) and *Cn( r )/ cr L ^ ie measurements in 30 /r _1 Mpc < r < 
78 h~ l Mpc. We choose this fitting range because there are less galaxy pairs at larger 
separations due to the subvolume size, and the nonlinear effect becomes too large for our 
SPT predictions to be applicable at smaller separations. For the two-point function, we 
take the Fourier transform of m 

Pg(k) = b\[P t (k) e - iV - + A MC P uc (k)\ , (5.13) 


where b\ is the linear bias, Pi(k) is the linear power spectrum, Amc is the mode coupling 
constant, and 


Puc(k) = 2f fS( 9 )fj(|k - q|)[F 2 (q.k - q)] 2 . 


(5.14) 
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5. Measurement of position-dependent correlation function 




b 







Figure 5.2: (Top left) £(r) of the mocks in real space. The gray lines show individual 
realizations, while the dashed line shows the mean. The black solid line shows the best¬ 
fitting model. (Top right) y 2 -histogram of the 600 mocks jointly fitting the models to £(r) 
and i^n(r )/ o\ in real space. The dashed line shows the ^-distribution with d.o.f.=36. 
(Bottom left) iCi( r )/ a i °f the mocks in real space for 220 /r _1 Mpc subvolumes. (Bottom 
right) Same as the bottom left panel, but for 120 h~ 1 Mpc subvolumes. 
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where iCi,SPT(^) and *Ci,fe 2 ( r ) are computed from eq. (5.3) with eq. (2.40) and eq. (2.42), 
respectively, and a\ l is the variance of the linear power spectrum computed from eq. (2.43), 
using the subvolume sizes of L = 220 and 120 h~ l Mpc and the redshift of z = 0.57. Note 
that the size of the subvolumes affects the values of a'j t . We determine L by first measuring 
b\ using the real-space two-point function of the entire survey, and then find L such that 
= o 2 l assuming the cubic top-hat window functiorj^J We find that these values 
(L = 220 and 120 h -1 Mpc) agree well with the cubic root of the total survey volume 
divided by the number of subvolumes, to within a few percent. 

We fit the models to £(r) and *Cn( r )/ c7 I of both subvolumes simultaneously by mini¬ 
mizing 

(5.17) 


x 2 = 




where C -1 is the inverse covariance matrix computed from the 600 mocks, Di and M t 
are the data and the model in the i th bin, respectively. The models contain three fitting 
parameters bi, b 2 , and A M c- 

The models computed with the mean of the best-fitting parameters of 600 mocks are 
shown as the black solid lines in figure 5.2 The best-fitting parameters are b\ = 1.971 ± 
0.076, b 2 = 0.58 ± 0.31, and A M c = 1-44 ± 0.93, where the error bars are 1-cr standard 
deviations. The agreement between the models and the mocks is good, with a difference 
much smaller than the scatter among 600 mocks. Upon scrutinizing, the difference in £(r) 
is larger for larger separations because the fit is dominated by the small separations with 
smaller error bars. On the other hand, for i£(r)/er| the agreement is good for two sizes of 
subvolumes at all scales of interest. This indicates that the SPT calculation is sufficient to 
capture the three-point function of the mocks in real space. 

The data points in figure 522 are highly correlated. To quantify the quality of the fit, 
we compute the y 2 -histogram from 600 mocks, and compare it with the ^-distribution 
with the corresponding degrees of freedom (d.o.f.). There are 13 fitting points for each 
measurement (£(r) and two sizes of subvolumes for iCn( r )/ CT i) an d three fitting parameters, 
so d.o.f.=36. The top right panel of figure 522 shows the y 2 -histogram. The dashed line 
shows the ^-distribution with d.o.f.=36. The agreement is good, and we conclude that 
our models well describe both £(r) and °f the mocks in real space. 

Our b\ is in good agreement with the results presented in figure 16 of [62], whereas 
our b 2 is smaller than theirs, which is ~ 0.95, by 1.2cr. This may be due to the difference 
in the bispectrum models. While we restrict to the local bias model and the tree-level 
bispectrum, |62l includes a non-local tidal bias Em 121 EH! and uses more sophisticated 
bispectrum modeling using the effective F 2 kernel [64^ [63]. In appendix [Dj we show that 
using the effective F 2 kernel and the non-local tidal bias in the model increases the value 
of b 2 , but the changes are well within the 1 -er uncertainties. Also, the differences of the 
goodness of fit for various models are negligible. 

The fitting range as well as the shapes of the bispectrum may also affect the results: 
the integrated correlation function is sensitive only to the squeezed configurations, whereas 


3 In principle, the shape of the window function also affects a\ ; , but we ignore this small effect. 
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5. Measurement of position-dependent correlation function 





Figure 5.3: Same as figure 5.2 but in redshift space. 


[62] includes more equilateral and collapsed triangle configurations. Understanding this 
difference merits further investigations. 


5.1.4 Measurements in redshift space 


Figure 5.3 shows the measurements of £(r) (top left) and iCi( r )/ (J L (220 h Mpc in the 
bottom-left and 120 h~ l Mpc in the bottom-right panels) of the mocks in redshift space. 
The gray lines show individual realizations, while the dashed lines show the mean. Similar 


to the analysis in section 5.1.3 , we fit the models in redshift space to the measurements in 
30 /r _1 Mpc < r < 78 h~ 1 Mpc. In this section, we use General Relativity to compute the 
growth rate, f(z ) ~ f2 m (z) 0 ' 55 , which yields f{z = 0.57) = 0.751. We shall allow / to vary 
when interpreting the measurements in the actual data. 

Since there is no baryonic acoustic oscillation feature on the scales we are interested in, 




































5.1 Measurement of PTHalos mock catalogs 


71 



^ 0.8 
0.7 
- 0.6 

- 0.5 

- 0.4 

- 0.3 
| 0.2 



Figure 5.4: Correlation matrix estimated from 600 mocks in redshift space. The figure 
shows o\ and °f 220 h~ l Mpc subvolumes from bin 0 to 13, o\ and of 

120 h~ x Mpc subvolumes from bin 14 to 27, and £(r) from bin 28 to 40. 


we model the redshift-space two-point correlation function as 

£,gA r ) = b l + -^McCmcM] K , 


where £,i,a v {r) an d Cmc(*") are given in eq. (5.15) and 


K = l + \l3 + \/ 


(5.18) 


(5.19) 


is the Kaiser factor with (3 = f /b\ [80]. As we do not include the subdominant term 
proportional to b 2 in the two-point function, it only gives constraint on b \, which we 
can then use to break the degeneracy with 62 in the integrated three-point function. We 
find that this simple modeling yields unbiased b\ and fulfills the demand. We calculate the 
redshift-space integrated three-point function by first evaluating the integrated bispectrum 
using SPT at the tree level (the explicit expression of the redshift-space bispectrum is given 


in appendix A.2), compute its one-dimensional Fourier transform, as eq. (5.4), and then 


a 


L,l 


correct for the boundary effect. The a\ of the mocks in redshift space agrees with b\K 
to percent level. The redshift-space models thus contain, as before in real space, the three 
fitting parameters, b \, b 2 , and Amc- We then simultaneously fit £(r) and iQj^r) j of 


both subvolumes by minimizing eq. (5.17). Figure 5.4 shows the correlation matrix {C t] in 


y 2 , normalized by yjCuCjj) estimated from the 600 mocks in redshift space. Because we 
normalize the integrated three-point function by the covariance between iQiSr) j and 
cr| is negligible. On the other hand, the covariances between i^^r)/a 2 L and £(r), between 
£(r), and between iCi( r )/ cr i f° r t wo sizes of subvolumes are significant. 

The models computed with the mean of the best-fitting parameters of 600 mocks are 
shown as the thick solid lines in figure 5.3. The best-fitting parameters are b\ = 1.931 ± 
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5. Measurement of position-dependent correlation function 



aV g[<mock 

var 

2 

mock. 

2 

^L,data 

220 h- 1 Mpc 

4.6 x 10" 3 

5.6 x 10~ 4 5 

4.9 x HT 3 

120 h- 1 Mpc 

2.4 x 10~ 2 

1.3 x HT 3 

2.5 x 10~ 2 


Table 5.1: Measurements of cr| of the mock catalogs and the BOSS DR10 CMASS sample. 


0.077, b 2 = 0.54 ± 0.35, and A MC = 1.37 ± 0.82. The agreement between the models and 
the measurements in redshift space is as good as in real space. 

Again, our b\ is in good agreement with the results presented in figure 16 of [62], whereas 
our b 2 is smaller than theirs, which is ~ 0.75, but still well within the 1-a uncertainty. As 
noted in section [5.1.3[ the adopted models of the bispectrum are different. In appendix |D| 
we show that using the effective F 2 and G 2 kernels and the non-local tidal bias in the 
model increases the value of b 2 . However, the changes are within the uncertainties, and the 
goodness of the fit is similar for different models. Thus, in this paper we shall primarily 
use the SPT at the tree level with local bias for simpler interpretation of the three-point 
function, but also report the results for the extended models. 


5.2 Measurement of the BOSS DR10 CMASS sample 


We now present measurements of the position-dependent correlation function from the 
BOSS DR10 CMASS sample]^] in NGC. The detailed description of the observations can be 
found in [3] [BJ. Briefly, the sample contains 392,372 galaxies over 4,892 deg 2 in the redshift 
range of 0.43 < z < 0.7, which corresponds to the comoving volume of approximately 
2 hr 3 Gpc 3 . We also weight the galaxies by wboss to correct for the observational system- 
atics. We follow section 5.1.1 to divide the observations into subvolumes. However, the 
observations have their own set of random samples, which are different from the ones of the 
mocks (the random samples of the mocks have slightly higher n and different h(z)), so we 
adjust the redshift cuts to be z cut = 0.5108, 0.5717 and z cut = 0.48710, 0.52235, 0.55825, 
0.60435 for the two resolutions, respectively. The resulting properties of subvolumes of the 
observations and mocks are similar. 

The mocks are constructed to match the two-point function of the observed galaxies, 
but not for the three-point function. Hence there is no guarantee that the three-point 
function of mocks agrees with the observations. We can test this using our measurements. 

The measurements of £(r) and '<Cr ( r ) / a L from the observations are shown as the solid 
lines in figure 5.5 the measurements of a\ is summarized in table 5.1 The measurements 
are consistent visually with the mocks within the scatter of the mock^J and we shall 
quantify the goodness of fit using y 2 statistics later. 

To quantify statistical significance of the detection of iQiXr)/(J 2 L and the goodness of fit, 


4 Catalogs of galaxies and the random samples can be found in http://www.sdss3.org 

5 These measurements of iCt( r )/ cr £ are done for one effective redshift. We compare iCL{r)/ a L °f the 
observations and mocks in different redshift bins in appendix [EJ finding that the observations and mocks 
are consistent at all redshift bins to within the scatter of the mocks. 
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Figure 5.5: Measurements of the BOSS DR10 CMASS sample (black solid lines). The gray 
lines show individual mocks in redshift space and the dashed line shows the mean of mocks. 
(Top left) £(r), (Bottom left) i£ L (r)/cr\ for 220 /r _1 Mpc subvolumes, and (Bottom right) 
*Ci( r )/°'i f° r 120 h Mpc subvolumes. (Top right) y 2 -histogram of the 600 mocks jointly 
fitting the three amplitudes to £(r) and iCi( r )/ (J I i n redshift space. The dashed line shows 
the x 2 -distribution with d.o.f.=38. The y 2 value measured from the BOSS DR10 CMASS 
sample is 49.3. 


we use the mean of the mocks as the model (instead of the model based on perturbation 
theory used in section 5.1.4), and fit only the amplitudes of Kl(j) f £(r), and a\ to 
the observations and the 600 mocks by minimizing eq. (5.17). Specifically, we use 0*(r) = 
4iO t mock (r) as the model, where Oi(r) = iQiir)/ a 2 h , 0 2 (r ) = £(r), and 0 3 = a 2 , with the 
amplitudes Ai, A 2 , A 3 . 

Table 5.2 summarizes the fitted amplitudes. The 1-cr uncertainties and the correlations 
are estimated from the 600 mocks. Since we normalize iCi,( r ) by a L-> the correlation between 
Ai and A 3 is small. On the other hand, A 2 and A 3 are correlated significantly because < 7 2 
is an integral of the two-point function [eq. (2.43) 
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5. Measurement of position-dependent correlation function 



Ai A 2 A 3 

1 -cr error 
best-fit (DR10) 

0.12 0.03 0.04 

0.89 1.02 1.08 



(Ai,A 2 ) (Ai,A 3 ) (A 2 ,A 3 ) 

corr 

0.34 0.09 0.36 


Table 5.2: Results of fitting the amplitudes: A\ is iCn( r )/ cr L A 2 is £(r), and A 3 is cr|. 
(Left) The 1-cr uncertainties of the amplitudes estimated from the mocks, and the best¬ 
fitting amplitudes of BOSS DR10 CMASS sample with respect to the mean of the mocks. 
(Right) The correlation coefficients of the amplitudes. 


Comparing the BOSS DR10 CMASS sample to the mean of the mocks, we find that 
i( L (r)/al is 1 -cr lower, £(r) is unbiased (by construction of the mocks), and cr 2 is 2 -cr higher. 
The result of A\ for the data is driven by the correlation between different separations of 
iQi^r) f o\. On the other hand, the result of A 3 is driven by the positive correlation between 
£(r) and o\. While o\ of the data for two subvolumes are larger than that of the mocks 
but still at the boundary of the variances (see table 5.1), it requires an even higher A 3 to 
minimize y 2 when we jointly fit the three amplitudes. The fact that A 3 is larger than A 2 
is also possibly due the contributions to a 2 from small separations (including stochasticity 
at zero separations), where the mocks were not optimized. We find A± = 0.89 ± 0 . 12 , i.e., 
a 7.4cr detection of the integrated three-point function of the BOSS DR10 CMASS sample. 

In order to assess the goodness of fit, we use the distribution of y 2 , a histogram of 
which is shown in the top right panel of figure 5.5 In total there are 41 fitting points (13 
fitting points for £(r) and two sizes of subvolumes for KL(r)/(? 2 L , and two fitting points 
for a 2 L ) with three fitting parameters, so d.o.f.=38. The y 2 value of the observations is 
49.3, and the probability to exceed this y 2 value is more than 10%. Given the fact that 
the mocks are constructed to match only the two-point function of the observations, this 
level of agreement for both the two-point and integrated three-point correlation functions 
is satisfactory. 


5.3 Interpretation for the measurement of the inte¬ 
grated three-point function 


What can we learn from the measured i^(r)/cr|? In section 5.1.4, we show that the pre¬ 
diction for i( L (r)/a1 based on SPT at the tree-level in redshift space provides an adequate 
fit to the mocks to within the scatter of the mocks; thus, we can use this prediction to 
infer cosmology from i<%(r)/cr 2 • Note that any unmodeled effects in the integrated three- 
point function such as nonlinearities of the matter density, nonlocal bias parameters, and 
redshift-space distortions beyond the Kaiser factor, will tend to bias our measurement of 
cosmological parameters based on i( L (r). We will discuss caveats at the end of this section. 

Since the linear two-point and the tree-level three-point functions are proportional to 
erf and < 7 g, respectively, and a 2 is proportional to erf, the scaling of the redshift-space 
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baseline 

eff kernel 

tidal bias 

both 


0.41 ±0.41 

0.51 ±0.41 

0.48 ±0.41 

0.60 ±0.41 


Table 5.3: Best-fitting b 2 and their uncertainties for BOSS DR10 CMASS sample for the 
extended models. The detailed description of the extended models is in appendix |D} 
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(5.20) 


where “fid” denotes the quantities computed with the fiducial value of cr 8 . Note that 
^Mc( r ) is proportional to <7g because it is an integral of two linear power spectra (see 
eq. (5.14)). Since £z, CT „(r) dominates the signal, the parameter combinations bias and 


K = 1 + 2/3/3 + (3 2 /5 are degenerate in the two-point function. That is, the amplitude of 
the two-point function measures only (bias) 2 ± §(bias)(fas) + |(/cr 8 ) 2 . This degeneracy 
can be lifted by including the quadrupole of the two-point function in redshift space. See 
[ 150 , 116411151 i HO] for the latest measurements using the BOSS DR11 sample. 


As for the three-point function, figure 2.2 shows that the bf and b\b 2 terms are com¬ 
parable for b\ b 2 . This means that, at the three-point function level, the nonlinear 
bias appears in the leading order, so the amplitude of the three-point function measures 
a linear combination of b\ and b 2 . This provides a wonderful opportunity to determine 
b 2 . The challenge is to break the degeneracy between b 2l bi, /, and a 8 . For this purpose, 
we combine our results with the two-point function in redshift space and the weak lensing 
measurements of BOSS galaxies. We take the constraints on bia 8 (z = 0.57) = 1.29 ± 0.03 
and f(z = 0.57 )a 8 (z = 0.57) = 0.441 ± 0.043 from table 2 of |130j . To further break 
the degeneracy between fq, /, and a 8 , we take the constraint on a 8 = 0.785 ± 0.044 from 
[ 1147117] . where they jointly analyze the clustering and the galaxy-galaxy lensing using the 
BOSS DR11 CMASS sample and the shape catalog from Canada France Hawaii Telescope 
Legacy Survey. 

We assume Gaussian priors on bia 8 , fa 8 , and a 8 with the known covariance between 
bia 8 and fa 8 . The cross-correlation coefficient between bia 8 and fa 8 is —0.59, as shown in 
figure 6 of [ 130] , We then run the Markov Chain Monte Carlo with the Metropolis-Hastings 
algorithm to fit he model eq. (5.20) to the observed /<C(r)/cr|. We find b 2 = 0.41 ± 0.41, 


and the results for the extended models are summarized in table 15.31 

The value of b 2 we find is lower than the mean of the mocks, 6™ ock = 0.54 ± 0.35. The 
difference is mainly due to two reasons. First, the amplitude of the integrated three-point 
function of the observations is lower than that of the mocks by 10% (Ai = 0.89 ± 0.12). 
Second, the priors from the correlation function and lensing constraint bi to be close to 
2.18, which is larger than that of the mocks, 5™ ock = 1.93. Thus, it requires a smaller b 2 
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5. Measurement of position-dependent correlation function 


to fit the three-point function. The argument is similar for the extended models. Note, 
however, that the nonlinear bias of the data is still statistically consistent with the mocks. 

Let us conclude this section by listing three caveats regarding our cosmological inter¬ 
pretation of the measured integrated three-point function. 


1. The models we use, eq. (5.20), are based on tree-level perturbation theory, the lowest 
order redshift-space distortion treatment, as well as on the local bias parametrization. 


While this simple model describes the mocks well, as shown in section 5.1.3 and 5.1.4 
we discuss in appendix [D] that using the effective F 2 and G 2 and the non-local tidal 
bias brings b 2 closer to that of |62j. We, however, find similar goodness of fit for 
various models, and thus we cannot distinguish between these models. 


2 . Covariances between the integrated three-point function, monopole and quadrupole 
two-point function, and weak lensing signals are ignored. This can and should be 
improved by performing a joint fit to all the observables. 


3. The cosmology is fixed throughout the analysis, except for / and erg- In principle, 
marginalizing over the cosmological parameters is necessary to obtain self-consistent 
results, although the normalized integrated three-point function is not sensitive to 
cosmological parameters such as Vt m as shown in figure 4.6 


These caveats need to be addressed in the future work. 


5.4 Discussion and conclusion 

In this chapter, we have reported on the first measurement of the three-point function 
with the position-dependent correlation function from the SDSS-III BOSS DR10 CMASS 
sample. The correlation between the position-dependent correlation function measured 
within subvolumes and the mean overdensities of those subvolumes is robustly detected at 
7.4a. 

Both the position-dependent correlation function and the mean overdensity are easier 
to measure than the three-point function. The computational expense for the two-point 
function is much cheaper than the three-point function estimator using the triplet-counting 
method. In addition, for a fixed size of the subvolume, the integrated three-point function 
depends only on one variable (i.e., separation), unlike the full three-point function which 
depends on three separations. This property allows for a useful compression of information 
in the three-point function in the squeezed configurations, and makes physical sense because 
the integrated three-point function measures how the small-scale two-point function, which 
depends only on the separation, responds to a long-wavelength fluctuation, as introduced 
in chapter |4j As there are only a small number of measurement bins, the covariance matrix 
of the integrated three-point function is easier to estimate than that of the full three-point 
function from a realistic number of mocks. We have demonstrated this advantage in this 
chapter. 
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Of course, since this technique measures the three-point function with one long-wavelength 
mode (mean overdensity in the subvolumes) and two relatively small-wavelength modes 
(position-dependent correlation function), it is not very sensitive to the three-point func¬ 
tion of other configurations, which were explored in [62]. 

We have used the mock galaxy catalogs, which are constructed to match the two- 
point function of the SDSS-III BOSS DR10 CMASS sample in redshift space, to validate 
our method and theoretical model. We show that in both real and redshift space, the 
integrated three-point function of the mocks can be well described by the tree-level SPT 
model. However, the nonlinear bias which we obtain from the mocks is higher than that 
reported in [62], This is possibly due to the differences in the scales and configurations of 


the three-point function used for the analyses. As discussed in section A3, any unmodeled 
nonlinear effects in the redshift-space integrated three-point function of CMASS galaxies 
will tend to bias b 2 , and will bias this parameter differently than the measurement of [62] . 

Taking the mean of the mocks as the model, and treating the amplitudes of two- and 
three-point functions as free parameters, we find the best-fit amplitudes of K{r)/cr 2 L , £(r), 
and (j\ of the CMASS sample. With respect to the mean of the mocks, the observations 
show a somewhat smaller iC( r )/ cr L (^1 = 0-89 ± 0.12) and larger erf,, while the ensem¬ 
ble two-point function £(r) matches the mocks. Given that the mocks are generated to 
match specifically the two-point function of the BOSS DR10 CMASS sample within a cer¬ 
tain range of separations, the level of agreement between the observations and mocks is 
satisfactory. 

Finally, by combining the integrated three-point function and the constraints from the 
anisotropic clustering (bias and fag in [ 130] ) and from the weak lensing measurements (a$ 
in [117] ). we break the degeneracy between bi, b 2 , /, and a 8 - We find b 2 = 0.41±0.41 for the 
BOSS DR10 CMASS sample. The caveat of this result is that our model, eq. (5.20), relies 
on a rather simple model in redshift space as well as on the local bias parametrization. We 
leave the extension of the model to improved bias and redshift-space distortion modeling 
(especially in light of the comparison with the results in [62]) for future work. 

In summary, we have demonstrated that the integrated three-point function is a new 
observable which can be measured straightforwardly from galaxy surveys using basically 
the existing and routinely applied machinery to compute the two-point function, and has 
the potential to yield a useful constraint on the quadratic nonlinear bias parameter. More¬ 
over, since the integrated three-point function is most sensitive to the bispectrum in the 
squeezed configurations, it is sensitive to primordial non-Gaussianity of the local type (see 
section 2.1.4 and section 4.5 for more details), thereby offering a probe of the physics 
of inflation. We plan to extend this work to search for the signature of primordial non- 
Gaussianity in the full BOSS galaxy sample. 
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5. Measurement of position-dependent correlation function 



Chapter 6 

The angle-averaged squeezed limit of 
nonlinear matter fV-point functions 
and separate universe simulations 


In the previous chapters, we have introduced the position-dependent two-point statistics 
and shown the measurements from iV-body simulations and observations. In this chapter, 
we shall generalize the study to the angle-averaged squeezed limit of nonlinear matter 
iV-point functions, and demonstrate how to use the “separate universe simulations” to 
address this issue. Specifically, we consider the case of the squeezed limit such that there 
is a hierarchy between two large wavenumbers k and k', and N — 2 small wavenumbers 
k l5 • • • , k n . The configuration is sketched in figure |6.1 


The squeezed limit of matter iV-point functions has been the subject of a large body 
of work in the context of the so-called “consistency relations” for the large-scale structure 
[851 EH [831 m [1651 EH EH [371 EUl m G31 H2U. The contributions to iV-point 
functions in the squeezed limit are ordered by the ratio of wavenumbers ki/k , which is 
assumed to be much less than one. The lowest order contributions, up oc (fcj/fc) -1 when 
the iV-point function is written in terms of the overdensity <5, are fixed by the requirement 
that a uniform potential perturbation as well as a uniform velocity (boost) do not lead to 
any locally observable effect on the density held, as demanded by the equivalence principle 
[S3., IT? 1 1166L M\ ■ They are also referred to as “kinematic contributions”. 

The next order contribution, oc ( ki/k)° , is the lowest order at which a physical cou¬ 
pling of long- and short-wavelength modes happens. More precisely, the contributions at 
this order correspond to the impact of a uniform long-wavelength density or tidal per¬ 
turbation. When considering equal-time iV-point functions and subhorizon perturbations 
ki aH , the kinematic contributions disappear, and the physical ( ki/k)° contributions 
are the leading contribution to the TV-point function in the squeezed limit. 


In this chapter, we disregard tidal fields, which leads us to first angle-average over the 
N — 2 small momenta (wavenumbers) in the iV-point function. Specifically, we consider 
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Figure 6.1: Sketch of the squeezed-limit configuration of N -point functions considered in 
this chapter. k 1; • • ■ , k„ denote the long-wavelength modes which are spherically averaged 
in eq. (6.1), while k and k' denote the small-scale modes which are allowed to be fully 
nonlinear. 


<Siv _2 defined through 

S N - 2 (k,k'\ h, ■ ■ ■ , k N - 2 ) = J ~ ^~ 2 (<^(k)5(k / )5(k 1 ) • • •5(k A r_ 2 )) / c (6.1) 


where ki are unit vectors and (5(ki) ■ ■ • <5(kAr))(, denotes the nonlinear connected matter Ap¬ 
point function with the momentum constraint (27r) 3 5D(ki -f • —h k^) dropped. Note that 
the momentum constraint fixes k' in terms of k and ki,..., k)v_ 2 . We now let k \,..., k^~ 2 
go to zero, and normalize the result by the nonlinear power spectrum P(k ) and the linear 
power spectra Pi(ki ) ■ ■ ■ Pi(k N _ 2 ) to obtain a dimensionless quantity: 


Rn- 2 (h) 


lim 

ki —^0 


<SN-2(k, k!\ ki, ■ ■ ■ , k N - 2 ) 
P( k)P l (k 1 )---P l (k N _ 2 ) 


( 6 . 2 ) 


Note that in this limit, spatial homogeneity enforces k' = —k + 0(ki/k), so that (for 
statistically isotropic initial conditions) the right-hand-sight of eq. (6.2) depends only on 
k. 


I 11 appendix [F] (see also M), we show that R n (k ) correspond exactly to the power 
spectrum response functions, which quantify the change in the nonlinear matter power 
spectrum to an infinite-wavelength density perturbation. These response functions are 
defined as the coefficients of the expansion of the power spectrum in the linearly extrapolated 
initial overdensity Sl 0 : 


P(k,t\5 LQ ) = Y J -R n {k,t) \s L 0 D(t)] n P(k,t) , 
' n\ L J 


(6.3) 


n=0 


where P(k,t\Si,o) is the nonlinear matter power spectrum at time t in the presence of 
a homogeneous (infinite-wavelength) density perturbation, and D(t) is the linear growth 
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factor normalized to unity today. We have set Ro(k,t ) = 1 by definition. Thus, by 
measuring R nj we measure the angle-averaged squeezed limit (eq. (6.2)) of the nonlinear 
matter (n + 2)-point function. For n — 1, the response Ri describes the angle-averaged 
squeezed-limit bispectrum, as discussed in chapter [2] and chapter [4} 


Independently of the derivation of eqs. ( 6.2| ) (6.3), we also present accurate measure¬ 
ments of R n for n = 1, 2, 3 using IV-body simulations which do not rely on approximations 
Specifically, we resort to IV-body simulations with an external liomoge- 


in section 6.2 


neous overdensity imposed via the separate universe approach described in chapter [3] A 
flat FLRW universe with a homogeneous overdensity is exactly equivalent to a different, 
curved FLRW universe, so that IV-body simulations in this modified cosmology provide, 
in principle, the exact result for the response functions R n {k). This in turn corresponds to 
the exact (in the limit of infinite volume and resolution) measurement of the squeezed-limit 
IV-point function (eq. (6.2)). The measurements of R\{k) are presented in |97j. We shall 
extend to n = 2 and 3. 


Many semi-analytical approaches to nonlinear large-scale structure assume that non¬ 
linear matter statistics can be described as a unique function of the linear matter power 
spectrum, i.e. the power spectrum of initial fluctuations linearly extrapolated to a given 
time. In the context of consistency relations, this approximation has been studied in 
e.g. |1651 [M], This ansatz is motivated by the fact that in Einstein-de Sitter (flat matter- 
dominated universe), and to a very good approximation in ACDM, the perturbation theory 
predictions factorize into powers of the linear growth factor and convolutions of products 
of the initial matter power spectra and time-independent functions. Another way to phrase 
this ansatz is that nonlinear large-scale structure only depends on the normalization of the 
fluctuations at a given time, and not on the growth history. 


In the context of squeezed-limit IV-point functions, this ansatz can be tested quanti¬ 
tatively by comparing the outputs of separate universe simulations at a given time with 
simulations in which the initial amplitude of fluctuations is rescaled to match the linear 
power spectrum at the same time. The difference between these “rescaled initial ampli¬ 
tude” simulations and the separate universe simulations corresponds to the error made in 
the ansatz of assuming that the linear power spectrum at a given time uniquely describes 
nonlinear large-scale structure at the same time. For n = 1, it is studied in [98] and found 
that the two simulations differ in the nonlinear regime. A closely related test using the 
matter bispectrum is shown in [120] . We shall study this comparison in more detail and 


for n — 1,2 and 3 in section 6.4 


This chapter is organized as follows. In section [63| we develop semi-analytic predictions 

e describe the methodology of performing 
we compare the measurements from N- 


for the power spectrum response. In section 6.2, we describe the methodology of performing 
the separate universe simulations. In section [673 


body simulations to the semi-analytic predictions. In section 6.4, we compare the rescaled 


simulations to the separate universe simulations. We conclude in section 6.5 
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6.1 Power spectrum response 


We define the n-th order response function R n (k ) of the power spectrum as the n-th 
derivative of the power spectrum with respect to the linearly extrapolated (or Lagrangian) 


overdensity 8l, normalized by the power spectrum. The definition, consistent with eq. (6.3), 


is 


Rn(k,t ) = 


1 d n P(k,t\8 L ) 


P(k) d[S L (t )]’ 


(6.4) 


<5l=0 


where 8i,(t) = 8loD{1)- In the following, we will frequently suppress the time argument for 
clarity. Analogously, one can define the power spectrum response functions with respect to 
the fully evolved (or Eulerian) nonlinear overdensity 8 p . Since we can expand the nonlinear 
overdensity in powers of 8l with known coefficients via the spherical collapse (see chapter [3]), 
the n-th order Eulerian response function is given by a sum of R m with m < n. Motivated 
by the relation eq. (6.2), we mainly consider the Lagrangian response functions. In the 


remainder of this section, we develop semi-analytic models for the response functions based 
on the separate universe picture. 


6.1.1 Separate universe picture 


The idea of absorbing an infinitely-wavelength overdensity perturbation 8 P into the back¬ 
ground modified (curved) cosmology is extensively discussed in chapter |3j and the response 
of the power spectrum to the overdensity is also discussed in chapter [4j (but restricted to 
the linear order). In this section, we shall first summarize results, and generalize the study 
to higher order response, i.e. n > 1. 

Because of the overdensity, the expansion slows down, and the overdense region behaves 
as a universe with positive curvature. The scale factor in the modified cosmology can be 
written as 

a(t) = a(t ) [1 + 8 a (t)} , (6.5) 


where the quantities in the modified cosmology are denoted with a tilde. Since the physical 
distance is the same in both cosmologies, eq. (6.5) implies the change of the comving 
distance to be 

a(t) 


x = 


a(t) 


— [1 + 5 a (t)] 


-i 


( 6 . 6 ) 


Furthermore, due to mass conservation, the fractional difference in the scale factor is related 
to the overdensity 8 P by 

1 + 8 p (t ) = [1 + 8 a (t)] 3 . (6-7) 


Using the separate universe picture, we consider the matter power spectrum in this 
patch just as that of a region with no homogeneous overdensity but properly modified 
cosmology. The modification of the cosmology is such that the shape of the linear power 
spectrum is unchanged, since the ratio of photon, baryon, and cold dark matter densities 
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is unmodified; moreover, the transfer function parameters are unchanged: Q m h 2 = Q r: 
and Qbh 2 = klfii 2 - Thus, only the growth of structure is affected. 


h 2 


The power spectrum that enters in the response given by eq. (6.4) is defined with 


respect to the background density and comoving coordinates of the fiducial cosmology. 
Hence, the power spectrum calculated for the modified cosmology has to be mapped to 
that with respect to the background density and comoving coordinates of the fiducial 
cosmology. As discussed in chapter [4[ this mapping yields the “reference density” and 
“dilation” contributions to the response. These can be calculated exactly at any scale k 
to any given order given the nonlinear matter power spectrum in the fiducial cosmology. 
That is, we do not need to run IV-body simulations to calculate these effects. They are 
thus merely “projection effects”, unlike the effect of the modified cosmology on the growth 
of structure, which requires simulations in order to provide an accurate estimate. 

Let us denote the power spectrum for the modified cosmology as P(k). Then, the 
reference density effect simply rescales the power spectrum as, 


P(fc) re ' S n "t , ( l + y 2 F(/c ) 


( 6 . 8 ) 


where the argument of P (k) is not modified. The dilation effect due to the change in the 
coordinates given by eq. (6.6) implies k —» k = (1 + 5 a )k and changes the power spectrum 
by (see chapter [4] for detailed derivation) 


p{k) dil = on (i + s a yp([i + s a } k ) 


Putting the two together and using eq. (6.7) yields 


P(k) = (l+6 p )P([l+5 a ]k) , 


(6.9) 


( 6 . 10 ) 


where all quantities are evaluated at some fixed time t. Note that one prefactor of (1 + 5 P ) 
cancels, since the effect of the increased density is partially canceled by the corresponding 
decrease in physical volume. As derived in chapter [3j for an Einstein de-Sitter fiducial 
universe (and to high accuracy in ACDM) 5 a {t ) and 5 p {t ) have series solutions of the form 


OO 


3 a (f) y J 


3 L0 D(t) 


oo 


n =1 


S p (t) = ^f n 

n= 1 


Slo D(t) 


( 6 . 11 ) 


where D(t) = D{t)/D(tfi) is the fiducial growth factor normalized to one at the epoch to to 


which we extrapolate 8lo = 5c(h))- Note that the values of e n and f n are given in eq. (3.35) 


and eq. (3.36), respectively. 


The third contribution to R n comes from the effect of the modified cosmology on the 
growth of structure, which as mentioned above is the physical contribution which requires 
At-body simulations for an accurate measurement. We thus define a set of growth-only 
response functions G n (k ) which isolate the nontrivial effect of the long-wavelength pertur¬ 
bation on the growth of small-scale structure, 


G n (k) = 


1 d n P{k) 


P(k) dS r j 


( 6 . 12 ) 


S L =0 
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That is, G n are defined as R n without the contributions from the reference density and 


dilation given by eq. (6.10). This definition is an extension of the similar decomposition for 


n — 1 shown in [147 l ; 84l [T5l [98] . Thus, the formula for the power spectrum (with respect 
to global coordinates) in the presence of a long-wavelength overdensity is given by 


P(A:|<5 L ) = (l + <y 


1 +E m 


n=l 


(6.13) 


- k=(l+8 a )k 

Clearly, by the Leibniz rule, at any given order n the total or “full” response R n (k ), 


eq. (6.4), is composed of the functions G m (k ) and the numbers e m , f m with 1 < m < n, 
where the e m multiply derivatives of Gi(k) and P(k) with respect to k (up to the n th 
derivative). Specifically, the first three full response functions are given by 


k P'(k) 

Ri{k)=h + ei —U + Gl {k) , 

P{k) 

R 2 (k) f , kP\k) , 2 k 2 P"(k) , G 2 (k) , , kP\k) 

-^- /2 + e2 W + ei W + ^ + /iei W 

kP'(k) 

+ fiGxik) + e 1 —^G 1 (k) + ei kG[{k) , 


(6.14) 


(6.15) 


Rs(k) kP'(k) , , , G 3 (k) , kP'(k) , f G 2 (k) , , kP\k) 

— JlGi(k)ei ^ h J3 H---h e 3 ~FTrT\ ^ A—^-f ./l e 2 


6 


P(k) 
2 k 2 P”{k) 


6 


P(k) 
kP\k) 


+ fieiZ 2FW + /2Gl(A:) + /2ei P( k ) 


e\k 


cm , jkP’ik)^,^ , ^ep"\k) 


+ e 


er 


2 11 P(k) 

kP'(k ) G 2 (fc) 


-fcG'Afc) + e?- 


+ (/iTi + ^2 )kG\(k) + e 
k 2 P"[k ) 


P(k) 
2 k 2 G'{{k) 


6P(k) 


+ 2eie 2 - 


2P{k) 


„, jW kP'(k) ,k 2 P"(k) 
+ Gi(fc) e 2 v ' 


P(fc) 2 w V p O) 2P(fc) 

where the primes denote derivatives with respect to /c. 


(6.16) 


6.1.2 Linear power spectrum predictions 


We now evaluate eq. (6.13) for the simplest case, i.e., the response of the linear matter 


power spectrum. In linear theory, the growth is scale-independent and given by the linear 
growth factor. Thus, the growth-only response functions are scale-independent and just 
described by the linear growth factor in the modified cosmology D(t), 


G 


linear 


1 d n (D 5 


D 2 dS n L 


(6.17) 


Sl= 0 


As for S a (t) and 5 p (t), D also has the perturbative expansion in powers of Sl as 


D(t) = D(t) 1 1 + £ g„ 


(6.18) 


n=1 
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Thus, for an Einstein-de Sitter fiducial universe (and to high accuracy in ACDM), the 
linear response functions are simply constants. Inserting the values of g n from eq. (3.49), 
we obtain 

r r iinean _ [ 26 3002 240272 197919160] 

*■ n ~ ’ 1323 ’ 43659 ’ 11918907 J 

eq. (6.13) evaluated for the linear matter power spectrum Pi(k,t ) then becomes 


(6.19) 


Pi(k,t\S L ) = [1 + 5 P (£)] 

Inserting the series expansions, we obtain 


m 

D(t) 


^,fid([l + S a (t)]k,t) . 


( 6 . 20 ) 


Pl(k,t\5 L0 ) = 1 + J2fn[SLoD(t)] n 1 + J^9n \SLoD(t) 


72—1 


72 — 1 


x Pi 


Lfid 


l + J2 e ^L 0 D(t)]' 


72—1 


M • 


( 6 . 21 ) 


Eq. (6.21) allows for a consistent expansion in 5lo■ Specihcally, d n Pi(k)/d5™ 0 is given by 
the ri-th order coefficient in this expansion, multiplied by n!. 


6.1.3 Nonlinear power spectrum predictions 

Beyond the linear matter power spectrum, the growth coefficients G n will become scale- 
dependent functions G n (k). Consider now what standard perturbation theory (SPT) pre¬ 
dicts. The power spectrum prediction is given by a series 

P SPT (k ) = Pi{k) + P 1 ~ loop (k) + P 2 ~ loop (k) + • • • , (6.22) 

where P n_lo °p scales as [Pi] n - In an Einstein-de Sitter universe, one can show (e.g., [17] ) 
that the time- and scale-dependence of each order in perturbation theory factorizes, so 
that one can write 


P SPT (k, t) = D 2 (t)P l (k, t 0 ) + D 4 (t)P 1 - lo °P(fc, t 0 ) + b Q P 2 ~ lo ° p (k, t 0 ) + 


(6.23) 


where P n loop (k,to) is a convolution of n factors of Pi(k,to) with time-independent co¬ 
efficients. While eq. (6.23) is only strictly correct in Einstein-de Sitter, it is used very 


commonly for ACDM as well, since departures from the exact result are typically of order 
1% or less, and since it simplifies the calculation significantly. Various variants of SPT, 
such as the renormalized perturbation theory (RPT) [4lj, share the same property. 


Eq. (6.23) allows for a very simple evaluation of the growth-only response: as discussed 


above, the shape of the linear power spectrum in the modified cosmology is uncha nged , 
and hence P SPT (k) can be simply evaluated by replacing the fiducial D(t) in eq. (6.23) 


with the modified one, eq. (6.18). This is equivalent to assuming that the entire late-time 
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cosmology dependence of the nonlinear matter power spectrum enters through the linear 


growth factor [1651 841 fT5lj . In section 6.4, we shall test this prescription to all orders in 
SPT calculations by performing simulations with a rescaled initial power spectrum. 

Apart from the SPT calculation, we can also apply this approximation to any prescrip¬ 
tion that maps a given linear power spectrum to a nonlinear one. In particular, we will 
show results for halofit [152] . In this case, where the dependence on the linear growth 
factor is not explicit, we instead compute the derivative with respect to the normalization 
of the linear power spectrum, 

A -4 A A . (6.24) 

dD dD da s 

which at the redshift considered yields the equivalent change of the linear matter power 
spectrum. This leads to 

D*AM -4 . (6.25) 

dD " s dd$ ' ' 

We use a five-point stencil with a step size of 0.75% in er 8 to compute numerically the 
derivatives with respect to Og. In conjunction with the change of the linear growth factor 
eq. (6.18), this allows us to compute the growth-only response G n (k ) for perturbation 


theory as well as fitting formulae of the nonlinear matter power spectrum. 


6.1.4 Halo model predictions 


In section 4.3.4 , we have derived the linear response of the power spectrum under the 
halo model framework, in which all matter is assumed to be contained within halos with 
a certain distribution of mass (given by the mass function) and density profile. In this 
section, we shall focus on generalizing the derivation to higher order responses for the halo 
model approach. 


by 


Adopting the notation in section 4.3.4 the halo model powe spectrum, P 


)HM 


is given 


P nM (k) = P 2h (k) + P lh (k) , P 2 \k) = [I\{k)] 2 P^k) , P lh (k ) = /°(M) , (6.26) 


where 


kM Y 


Imih, ■ ■ ■ , k m ) = j d 111 M 71 (In M) J b n (M)u(M\ki) • • • u(M\k m ) , 


(6.27) 


and n(lnM) is the mass function (comoving number density per interval in log mass), 
M is the halo mass, b n (M) is the n-th order local bias parameter, and u(M\k) is the 
dimensionless Fourier transform of the halo density profile, for which we use the NFW 


profile. One can find more details in section 4.3.4 


To derive how the power spectrum given in eq. (6.26) responds to a homogeneous 
(infinitely long-wavelength) density perturbation 5l, we consider the one-halo (-P lh ) and 
two-halo (P 2h ) terms separately. The key physical assumption we make is that halo profiles 
in physical coordinates are unchanged by 8l- That is, halos at a given mass M in the 
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presence of 5l have the same scale radius r s and scale density p(r s ) as in the fiducial 
cosmology. We will discuss this assumption later in this section. Given this assumption, the 
density perturbation 5l then mainly affects the linear power spectrum, which determines 
the halo-halo clustering (two-halo term), and the abundance of halos at a given mass. 

For the two-halo term, as the response of the linear power spectrum is derived in 
eq. (6.21), the remaining task is to consider the effect on I\. By the assumption that the 
density profile does not change in the presence of Sl, I\ only changes through the bias 
b\(M ) and the mass function n(lnM). The bias bm(M) quantifies the TV-th order response 
of the mass function n(ln M) to Sjj: 


b N (M ) = 


<9^77,(111 M) 


n(lnM) 


as? 


or 


d N n(ln M) 


d5 n l 


= b N (M)n(\n M) . 


(6.28) 


Thus, 

d N 


as? 




<5l=0 


M\ B N 

= I dhrM | — 1 [bi(M)n(\nM)\ * u(M\k) = Ii +1 (k) . (6.29) 


P / 

Note that in the large-scale limit, k —» 0, this vanishes for > 1 by way of the halo model 
consistency relation 


Ir(k) = 


d In M n(ln M) ( —^ b^{M) = 


1, N — 1 
0, N > 1 


For hnite k eq. (6.29) does not vanish, we thus have 


1 


ll(k,t\5 L0 ) = Y, -I? + \k,t)iD(t)6 L o] 


(6.30) 


(6.31) 


71—0 


Thus, the two-halo term in the presence of becomes 


P 2h (k, t\5 L0 ) = ll + Yfn[h 0 m] 


71=1 
OO ^ 


X 


\n =0 



e n [Sw D(t)Y 


(6.32) 


k, t 


Note that we recover the tree-level result given in eq. (6.21) in the large-scale limit, where 
only n — 0 of the third term in eq. (6.32) survives. Note also that in eq. (6.32) the 
dilation effect only enters in the linear , not 2-halo, power spectrum. This is a consequence 
of our assumption that halo profiles do not change due to the long-wavelength density 
perturbation. 

For the one-halo term, due to our assumption about density profiles, the only effect is 
the change in the mass function, which through eq. (6.28) becomes 

d N 


i"(k,k) = i?(k,k), 


as? 


(6.33) 






















6. The angle-averaged squeezed limit of nonlinear matter iV-point functions 
88 and separate universe simulations 


JLOl 


and thus 

P lh (k,t\5 L0 ) = T-J^k,k,t)[D(t)6 L 

n\ 

n =0 

Putting one-halo and two-halo terms together, we obtain 

p hm (mim = (i + ^/„Mwr) ( 1 + [< h 0 D(t ) 

n= 1 

\ 2 / r 

1 ....~ ' i+j2 e n[h 0 mr 


(6.34) 


n= 1 
oo 


X 




Z,fid 


vn=0 
oo 


k, t 

/ \ n= 1 

+ E^! / 2(^A:,t)[D(t)(5 i o] n . (6-35) 

n =0 

The contribution oc I™ +1 (for n > 0) is numerically much smaller than the other terms (see 
also the discussion in section 4.3.4). Since it is much smaller than the overall accuracy of 
the halo model description, we will neglect it in the following. This yields 

/ OO \ / OO 

P UM (k,t\5 L0 ) = 1 + Y J fn[tLvD(t)] n \ 1 + Y j9 nU L oD(t) 


n =1 


72—1 


x (Il{k,t)) 2 P lM 


1 + J2 e n[6LoD(t)Y 


72—1 


k, t 


j^^(k,k,t)ims w ] n ■ 


(6.36) 


n=0 


Explicitly, the Erst and second order full response functions are given by 

d\nPi(k,t) 


Rf M (k) = 
R™(k) = 


f 1 + 2g 1 +e 1 - 

Cf ffJ. IV 

2/2 + 2 /i^i + (/1 + 2 gi)ei 


P 2h (k,t) + I^(k,k,t) 

+ 2 g 2 + 4^2 


d In Pi(k, t) | o 2 


din k 


din Pi(k, t) 2 1 d 2 Pi(k,t ) 
2 dIn k 1 P d(lnfc) 2 


P 2h (M) + / 2 (^M) 


(6.37) 


We also derive the growth-only response functions in the halo model approach. Since 
the halo profiles are assumed fixed in physical coordinates, this means that we need to 


rescale the halo model terms, accordingly. Following our discussion in section 6.1.1 


we have k — (1 + S a )k, where k is the comoving wavenumber with respect to the modified 
cosmology. We then obtain 


growth only 


&,■■■ ,~k m ) = r„ 


k\ 


kr i 


physical \1 + S a (t) 1 + 5 a (t) 


(6.38) 
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Inserting this into eq. (6.35) and performing a series expansion of 5 a in 5l then allows us to 
derive the growth-only response functions G^ M (fc). Note that the NFW profile we assume 
is uniquely determined by the scale radius r s (M) for a halo of mass M, which enters the 


coefficients defined in eq. (6.27) in the combination kr s (M). Thus, it is easily possible to 
include a dependence of the scale radius r s (M), or equivalently the halo concentration, on 
the long-wavelength density in a similar way. We will leave this for future work. 


Quantitatively, the main contribution of the rescaling eq. (6.38) is from the one-halo 


term oc I£(k,k), i.e. the term in the last line of eq. (6.36). The rescaling of the other 
instances of only changes the response at the sub-percent level and we will neglect 
them in the following. We then obtain for the growth-only contribution to the halo model 
power spectrum 


P hm (M!M growt = only (l + 5>n S w D(t) ") [Il(k,t )] 2 P lfid (k,t) 


n= 1 


+ X] ^ J 2 [A(S W , t ) k, A(5 L0: t)k , t] [D(t)5 L0 ] n , (6.39) 


n =0 


where 


-l 


A(S L0 ,t)= l + X>n[<WW 


71=1 


(6.40) 


6.2 Separate universe simulations 

To test our semi-analytical models of the power spectrum response to the homogeneous 
overdensity, particularly for the growth of structure due to the change of the cosmology, 
we run separate universe simulations and measure the power spectrum response func¬ 
tions directly. In this section, we describe the details for performing the separate universe 
simulations. We shall first introduce the straightforward modifications, cosmological pa¬ 
rameters and the initial conditions, and then the non-trivial choices, comoving or physical 
coordinates. 

For usual iV-body codes, one needs to specify the cosmological parameters at present 
time and the output time t out normally specified by the scale factor. As discussed in 
chapter [ 3 J in the presence of the overdensity 5lo, the cosmological parameters at a(t 0 ) = 1 
are modified as 

Ho = Hq(1 + 5h) , &m — ^m(l + $h) 2 , = ^a(1 + 8h) 2 

/ ~ \ 1/2 

-i. wrlm SL °- (6 - 41) 

For the output time, because we want to compare the simulations at the same physical 
time, we need to determine the corresponding scale factor in the modified cosmology such 
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that a(t ou t) = [1 + ^a(tout)]- We can numerically solve the ordinary differential equation of 
5 a in the modified cosmologies eq. (3.18); alternatively we can numerically evaluate a(t out ) 
by 


'■a(iout) 


tout. 


da 


n a(t. out) 


da 


r 0,(tout ) [l+^a (tout)] 


da 


aH(a) 


aH(a) 


aH(a) 


(6.42) 


In order to generate the initial conditions for iV-body simulations of the modified cos¬ 
mologies, we need the linear power spectrum at the initial redshift. The initial power 
spectrum has to be generated for the cosmology [fi m , 12 a, H 0 \ with the same amplitude of 
the primordial scalar curvature perturbation A s as for the fiducial cosmology. Since the 
transfer function only involves the physical matter and radiation densities quantified by 
O m ih ( j and so on, it is identical in the modified and fiducial cosmologies. Therefore, the 
linear power spectra differ only through the difference in the linear growth factor. We use 
CAMB [96) [95J to compute the power spectrum of the fiducial cosmology at z = 0, and 
rescale it by 


D(di) 


D(d 

D{a 


1 ) 

1 ) 


(6.43) 


where a* is the scale factor for which the initial conditions are generated]]] Next, we 
generate a Gaussian realization of the density field following the initial power spectrum. 
The positions and velocities of the particles are computed by the second-order Lagrangian 
perturbation theory m- 

Given a fixed box size for the fiducial cosmology, there are two reasonable choices for 
the box sizes of the modified cosmologies. Either we match the respective comoving box 
sizes, i.e. the box size is 500 h/h in units of h Mpc comoving, or we choose the box 
sizes such that their physical sizes coincide with that of the fiducial simulation at one 
specific output time t out , i.e. 500 ha(t oni )/[hd(t oni )] in units of h~ l Mpc comoving. The 
former choice is adequate if we are interested in the power spectrum response functions 
at the same comoving wavenumber, i.e. without the “dilation” effect. By using the mean 
density of the separate universe cosmology as the reference density when computing the 
power spectrum, we are further removing the “reference density” effect and are left with 
the growth-only response. The results of these simulations are presented in section |6.3.1 


In order to measure the full response functions, we run simulations for which we match 
the physical box size. We focus on two different output times i out corresponding to £ = 0 
and z — 2 in the fiducial cosmology. As the physical size can only be matched at one 
specific time, we have to run a new set of simulations for each output time. The results of 


these simulations are presented in section 6.3.2 


Tn order to recover the correct linear power spectrum at low redshifts, we compute the growth functions 
(D and D) without taking radiation into account. This is because IV-body codes do not include the effect 
of cosmological radiation. In our procedure, we also neglect the effect of curvature on the transfer function 
at very low wavenumbers k ~ \/\K\, since terms of similar order are neglected in the Poisson equation 
used in IV-body codes. For sub-horizon box sizes these effects are completely negligible. 
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Finally, we summarize the common features for the separate universe simulations. 
All simulations are gravity-only simulations and are carried out using the Tree-PM code 
Gadget-2 |153j . The starting redshift is z = 49 (a, = 0.02), and the particle load for each 
simulation is 512 3 . For the fiducial cosmology (£ i0 = 0), we choose a flat ACDM cosmol¬ 
ogy with cosmological parameters consistent with the current observational constraints: 
VL m = 0.27, h = 0.7, f \h 2 = 0.023, n s = 0.95, as = 0.8, and a comoving box size of 
500 h Mpc -1 . 

We simulate separate universes corresponding to the linearly-evolved present-day over¬ 
densities of 5lo = 0, ±0.01, ±0.02, ±0.05, ±0.07, ±0.1, ±0.2, ±0.5, ±0.7, and ±1. Then, 
for the separate universes, the Hubble constant and the curvature fraction vary between h: 
0.447 to 0.883 and f Ik- —2.45 to 0.372, respectively. The physical densities Vt m h 2 ^ Cl\h 2 , 
and f \h 2 as well as n s and the amplitude of the primordial curvature power spectrum 
remain the same. 


6.3 Results of separate universe simulations 

For the power spectrum computation, we first estimate the density contrast S(x) on a 
1024 3 grid using the cloud-in-cell mass assignment scheme, then apply a Fast Fourier 
transform, and angular average the squared amplitude |<5k | J ■ The density contrast 4(x) = 
p(x)/p — 1 describes the overdensity with respect to the reference density p. When we 
are interested in the growth-only response function, p is equal to the mean density of the 
separate universe. When we compute the full response function, p is equal to the mean 
density of the fiducial cosmology. Similarly, for the growth-only response, distances are 
measured using the comoving coordinates of the respective cosmologjj^J whereas, for the full 
response, the power spectrum is always measured in comoving coordinates of the fiducial 
cosmology. 

We only report results up to a maximum wavenumber of 2 h ^ 1 Mpc. A convergence 
study with simulations with 8 times lower mass resolution shows differences in G i, G 2 and 
Gs of only 1 (3) to 5 (10) percent at z = 0 (z = 2) up to that wavenumber, where the 
deviations increase from the linear response function to the higher-order response functions. 
The results for the full response functions Ri, i? 2 and R 3 are converged to an even better 
degree. We therefore expect that the simulation results presented here converge to a sub¬ 
percent to a few percent level. 

I 11 order to compute the first three response functions, we fit a polynomial in 5l to the 
fractional difference in the measured power spectrum A k($L) = P{k\5 l)/P{ k\5 l = 0) — 1 
for each k- bin. For the fit, we only include results from separate universe simulations 
with d/-,(Gut)| A 0.5 and use a polynomial with degree 6 to be unbiased from higher- 
order response functions. As the random realization of the initial density held is the same 
across different 5l values, the corresponding power spectra are strongly correlated. By 
considering the ratio, or the relative difference, of two power spectra a large fraction of 
the noise cancels. However, for the same realization, the measured fractional differences 


2 Note, however, that the unit of length is always h 1 Mpc, where h corresponds to the fiducial cosmology. 
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A k{$L) are still correlated over 8 l■ As the number of realizations (16) is not large enough 
to reliably estimate the covariance between different 8 l values, we cannot include this 
correlation in the polynomial fitting. Instead, we construct quasi-decorrelated samples 
of A k{8i) by randomly choosing a realization for each 8 l value. Fitting many of those 
subsamples allows for a robust error estimation of the derived response functions. 


6.3.1 Growth-only response functions 


Figure 6.2 shows the first three growth-only response functions measured from the simu¬ 
lations at z = 0 (left column) and z — 2 (right column). These correspond to the fully 
nonlinear squeezed limit bispectrum (three-point function), trispectrum (four-point func¬ 
tion) and five-point function. The small wiggles in the growth-only response functions 
result from the damping of the baryon acoustic oscillations (BAO), which depends on the 
amplitude of density fluctuations and thus on 8 l- 

Let us compare the simulation results to the theoretical predictions discussed in sec¬ 
tion 6.1 On sufficiently large scales, the perturbation theory predictions are the most 
accurate, as expected. At high redshift, the 1 -loop predictions best describe the results 
overall. The 1-loop predictions also show a BAO damping effect. At z = 0, the growth-only 
response is captured best by the haloht prescription (in case of G\) or the halo model (in 
case of G 2 , G 3 ). We see that the haloht prescription describes the simulation results of the 
linear response well at both redshifts, but performs significantly worse for the higher-order 
response functions. The BAO damping effect is essentially absent in both haloht and halo 
model predictions. Overall, none of the models is able to accurately describe the simulation 
data in the nonlinear regime, with discrepancies at z = 0 ranging from 20 % in the best 
case to a factor of several. These discrepancies are not surprising given that we are looking 
at scales beyond the validity of perturbation theory and at higher Appoint functions for 
which the semi-analytical approaches were not tuned. 


The halo model prediction does not asymptote exactly to the linear result in the k —> 0 
limit. This is because the one-halo term asymptotes to a white noise contribution in this 
limit, and since the one-halo term contributes to G n due to the dependence of the halo 
mass function on 8 l (section 6.1.4), this induces a correction to the linear prediction which 
contributes on large scales. Physically, this occurs because the halo model does not enforce 
momentum conservation of the matter density held. This issue can be fixed by introducing 
a “mass compensation scale” mg. 


The halo model predictions can be tuned to better match the simulation results by 
allowing for a dependence of the halo profiles on the long-wavelength density, which is 
expected on physical grounds (see also [98]). Specifically, if the scale radius of halos at 
fixed mass increases in the presence of a long-wavelength density perturbation, this lowers 
the peak in the response and thus could lead to better agreement with the simulations 
results. 
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Figure 6.2: The first three growth-only response functions of the power spectrum measured 
from the separate universe simulations at z = 0 (left) and z = 2 (right). The error bars 
show the statistical error derived by random resampling of the data (see text). For data 
points apparently without error bars, the statistical error is smaller than the size of a dot. 


6.3.2 Full response functions 


We now turn to the results for the full response functions, i.e. including the “dilation” and 
“reference density” effects. The results of the simulations and the model predictions are 
shown in figure [673 The oscillations in the response functions can be traced back to the 
BAOs in the power spectrum. The BAOs propagate to the response functions primarily by 























6. The angle-averaged squeezed limit of nonlinear matter iV-point functions 
94 and separate universe simulations 



Figure 6.3: The first three full response functions of the power spectrum measured from 
the separate universe simulations at z = 0 (left) and z = 2 (right). 


the “dilation” 
eqs. (6.14)— (6.16[)). 


effect, which yields derivatives of the power spectrum with respect to k (see 
The 1-loop perturbation theory predictions describe the simulation 
Mpc and k < 0.3 h~ l 


results accurately up to k < 0.15 h^ 1 Mpc and k < 0.3 /r _1 Mpc at z = 0 and z — 2, 
respectively. As the other theoretical models do not include the damping of the BAOs in 
the nonlinear power spectrum, they predict oscillations in the response functions which 
are too large. To improve the accuracy of those models around the BAO scale, one would 
need to put in the BAO damping by hand. In the nonlinear regime, none of the models is 
able to reproduce the simulation data. In principle, one could build a hybrid model for the 
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Figure 6.4: The first three Eulerian response functions of the power spectrum measured 
from the separate universe simulations (data points) at z = 0 (left) and z — 2 (right). The 
lines show the corresponding linear combinations of the Lagrangian response functions 
using the f n coefficients derived for the Einstcin-de Sitter universe (see eq. (6.44) and 
eq. g36h). 


full response by combining an accurate prediction of the nonlinear power spectrum of the 
fiducial cosmology and the growth-only response functions G n (k ) discussed in the previous 
section. However, we do not pursue this approach. 


6.3.3 Eulerian response functions 

So far, we have always considered the response to the linearly-extrapolated initial (La¬ 
grangian) overdensity 5^. We now consider the corresponding response to the evolved 
nonlinear (Eulerian) overdensity S p . Using the expansion derived for the Einstcin-de Sitter 
universe, eq. ( |3.36[ ), we find 

^Eulerian ^ = ? 

^Eu'enan^ = ^ _ 2 f 2 R^ k ) , 

R$ ule ™ Q (k) = R 3 (k ) - 6 f 2 R 2 (k) + 6 (2/| - / 3 ) R 1 (k) . (6.44) 

In figure |6~4| we compare the directly measured Eulerian response functions with the appro¬ 
priate linear combinations of the measured Lagrangian response functions. The agreement 
is excellent as expected, especially at high redshift at which the ACDM universe is very 
well approximated by the Einstein-de Sitter universe. 

Interestingly, the higher-order Eulerian response functions are much smaller than in 
the Lagrangian case. That is, the response of the nonlinear matter power spectrum to a 
uniform nonlinear final-time density S p is close to linear. This is most likely due to the fact 
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that the growth-only response functions are subdominant compared to the rescaling and 


reference density contributions, especially at higher order. In this case, eq. (6.10) implies 
a close to linear scaling with S p . 


6.4 Simulations with rescaled initial amplitudes 


All models for the growth-only response functions that we have presented in section 6.1 


and section 6.3.1 are based on the approximation that we can trade the effect of 8 l for an 
appropriate change to the linear growth factor (or equivalently, the linear power spectrum). 
But how well does this approximation work? 

To investigate how well the effect of a homogeneous overdensity on the growth of 
structure can be modeled by a change in the amplitude of the linear power spectrum, we 
additionally run a set of simulations for which we always assume the fiducial cosmology 
but vary the amplitude of the initial power spectrum. Specifically, for each 5lq value for 
which we simulate a separate universe, we also simulate the fiducial cosmology with the 
initial power spectrum amplitude multiplied by D(to) 2 / D(to) 2 , where D(to) is the linear 
growth factor in the corresponding separate universe cosmology. 

Using these “rescaled-amplitude simulations”, we can explicitly test the approximation 
that Sl effects the growth of structure only through the change in the linear growth factor 
on all scales including the nonlinear regime. 


6.4.1 Comparison to separate universe simulations 


In figure |6.5[ we show the growth-only response functions measured from two different 
sets of simulations. In case of G i, this comparison was also shown in figure 6 of [98], 
and our results agree with theirs]^] The rescaled-amplitude simulations all assume the 
fiducial cosmology but vary the amplitude of the linear power spectrum used to initialize 
the simulations so as to match the linear power spectrum in the modified cosmology at the 
given output times [eq. (6.18) and eq. (6.24)]. On linear scales, these simulations thus agree 
with the separate universe simulations by construction. As the simulations share the same 
random realization of the initial density, the sample variance (noise in the upper panels) 
gets vastly reduced when considering the difference of the measured response functions, 
A G n = G^ escaled — ^separate, This j s shown in the lower subpanels of figure 6.5, where we 
have divided A G n by the corresponding linear growth-only response, i.e. the prediction in 
the k —y 0 limit. 


The differences seen in figure [6~5] are caused by the different growth history, which is not 
captured by the rescaling of the initial amplitude. Following the discussion in section |6. 1.3 


the commonly used SPT approach factorizing the growth factor and scale dependence 
assumes at all orders that a long wavelength density perturbation enters exclusively through 
the modified linear growth. Thus, even when calculated to all orders , the best that this 


3 Note that the in [ T2'0] a different comparison is performed using the time derivative of the nonlinear 
power spectrum in simulations of the fiducial cosmology. 
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SPT calculation could do is to reproduce the rescaled amplitude result in figure |6.5[ which 
deviates from the actual response at z = 0 by 10% at k ~ 0.5 h~ l Mpc and 20% at 
fc ~ 1 /r 1 Mpc for Gi, and significantly worse for the higher-order response functions. At 
z = 2 on the other hand, the rescaled-amplitude G\ matches the separate universe response 
to better than 10% even beyond k = 1 h Mpc, and for G 2 , G$ performs significantly 
better as well. 

There are two possible explanations for these discrepancies in the SPT context. First, 
using the SPT kernels derived for an Einstein-de Sitter universe (which have time-independent 
coefficients), with the ACDM linear growth factor replacing the Einstein-de Sitter a(t), 
could become highly inaccurate for ACDM at higher orders. Note that the same issue 
exists for a fiducial flat Einstein-de Sitter universe, since for Sl % 0 the quantity D m // 2 
is no longer 1 (d(Q m // 2 )/ d5 l = —5/21 [15]; see also the discussion in [120] ). There is no 
indication of such a strong effect at low orders in perturbation theory, where this approx¬ 
imation typically performs to better than a percent im Furthermore, it is found in [15] 
that a cancellation in the curvature contribution to the growth integral suppresses this 
effect. Finally, it is shown in [98] that the growth-only response of the power spectrum to 
a change in the Hubble constant while keeping Q m h 2 fixed follows the separate universe 
response very closely (figure 6 there). If the much larger discrepancies between separate 
universe response and rescaled amplitude response were due to the cosmology dependence 
of the SPT kernels, one would not expect this to be the case. Nevertheless, we do not 
claim to be able to rigorously exclude this possibility. 

The other possibility, more likely in our opinion, is that the discrepancy between 
rescaled amplitude and full separate universe simulations is due to effective non-perfect 
fluid terms, such as pressure and anisotropic stress, in the dark matter fluid [T3]. The 
effective fluid properties depend on highly nonlinear small scales which are not described 
by the Euler-Poisson system. Their value can depend on the growth history (as well as 
the power spectrum shape) thus leading to a discrepancy between rescaled amplitude and 
separate universe simulations. Assuming this interpretation is correct, figure 6J3 explicitly 
shows the breakdown of SPT on nonlinear scales as effective pressure, anisotropic stress 
and sound speed need to be included. Separate universe simulations can be used to mea¬ 
sure the response of these effective terms to a long-wavelength overdensity, which is crucial 
when modeling (N > 2)-point functions. The results shown in figure 6.5 are analogous 
to what has been found for the mass function of halos which is a key ingredient in the 
halo model description of the nonlinear matter density field. The mass function shows 
departures from being a simple function of the linear matter power spectrum at the 5-10% 
level lIS'3. [ 20 ]. 

In an Einstein-de Sitter cosmology with scale-invariant initial power spectrum Pi(k) oc 
k n , there is only one characteristic spatial scale at any given time, which corresponds to 
the scale at which the density field becomes order 1 m- Let us denote this wavenumber 
as &nl(£)- Then, the response functions have to follow a universal function of k/k^Jt), i.e. 
G m (k,t ) = G m (k/ /cnl) (keeping the index of the initial power spectrum fixed). Thus, in 
this specific case, separate universe simulations and rescaled-amplitude simulations will give 
exactly the same result when compared at fixed k/k^L- The departures shown in figure [ti~5] 
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can thus be seen as a consequence of the ACDM background and the departure from scale- 
invariance of the initial power spectrum. It would be interesting to disentangle the two 
effects, e.g. by performing separate universe simulations in ACDM with scale-invariant 
initial conditions. We leave this for future work, but point out that when plotting the 
differences shown in the lower panel of figure 6.5 as a function of we still find a 

factor of several difference in the z = 0 and z — 2 results. 


6.5 Discussion and conclusion 

In this chapter, we described in detail the procedure for performing A r -body simulations 
with the separate universe technique. Using the separate universe simulations, we compute 
the response of the nonlinear matter power spectrum to a homogeneous overdensity super¬ 
imposed on a flat FLRW universe. The response functions we computed give the squeezed 
limits of the 3-, 4-, and 5-point functions, in which all but two wavenumbers are taken to 
be small and are angle-averaged. By virtue of the separate universe technique, we reach 
an unprecedented accuracy of these nonlinear matter Appoint functions. 

The response function consists of three parts: changing the reference density with 
respect to which the power spectrum is defined; rescaling of comoving coordinates; and the 
effect on the growth of structure. The former two effects can be calculated trivially, whereas 
the third one requires separate universe simulations. We have compared the simulation 
results with analytical and semi-analytical results, in particular standard perturbation 
theory (SPT), the empirical fitting function haloht, and the halo model, finding that SPT 
typically yields the best results at high redshifts. The fitting function and halo model, while 
qualitatively describe the trends seen in the response functions, give a poor quantitative 
description on nonlinear scales. 

A fundamental assumption of all of the analytical and semi-analytical methods used in 
this chapter, including standard perturbation theory at any order, is that nonlinear matter 
statistics at a given time are given solely by the linear power spectrum at the same time, 
and do not depend on the growth history otherwise. As was done in [98] for the response 
function for n — 1, we were able to test this assumption for n = 2 and 3 quantitatively 
by comparing the separate universe simulations with simulations with a rescaled initial 
power spectrum amplitude. We find that this assumption fails at the level of 10% at 
k ~ 0.2 — 0.5 h^ 1 Mpc for 5- to 3-point functions at z = 0. The failure occurs at higher 
wavenumbers at z = 2. In the context of SPT, this may signal a breakdown of the perfect 
fluid description of the dark matter density held at and beyond these wavenumbers. In 
other words, even if computed to all orders, SPT (and its variants such as RPT [41] ) 
fails to describe the nonlinear structure formation beyond these wavenumbers. Therefore, 
our results yields a quantitative estimate for the scales at which effective fluid corrections 
become important in the bispectrum and higher Appoint functions, and at which one 
should stop trusting pure SPT calculations. 

Finally, we point out that the approach presented here can be augmented to mea¬ 
sure more general squeezed-limit Appoint functions, by including the response to long- 
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wavelength tidal fields and by considering the response of small-scale n-point functions in 
addition to the small-scale power spectrum considered here. 
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Figure 6.5: Comparison of the growth-only response functions G i, G 2 , G 3 (top to bottom) 
measured at z = 0 (left column) and z = 2 (right column) from one realization of the 
separate universe simulations and from the same realization simulated by varying the 
initial amplitude. The bottom subpanels show the difference, A G n = C^ escaled — G® eparate , 
divided by the response of the linear matter power spectrum, G|) near . 












































Chapter 7 

Summary and outlook 


In this dissertation, we have proposed and developed in detail a new observable, position- 
dependent power spectrum, to extract the squeezed-limit bispectrum of the large-scale 
structure by measuring the correlation between the position-dependent two-point statis¬ 
tics and the mean overdensity in the subvolumes of a survey volume. Since this new tech¬ 
nique requires essentially measurements of the two-point statistics and mean overdensity, 
it sidesteps the complexity of the traditional three-point function estimation. 

The correlation between the position-dependent two-point statistics and the mean over¬ 
density can be regarded as how the small-scale structure formation responds to a long- 
wavelength mode. In chapter [3j we have shown that the long-wavelength overdensity 
compared to the scale of interest can be absorbed into the background cosmology, and 
the small-scale structure formation evolves in the corresponding modified cosmology. This 
separate universe approach thus provides an intuitive way to model the squeezed-limit, bis¬ 
pectrum, i.e. the response of the small-scale power spectrum to a long-wavelength mode. 

In chapter |4| we have measured the position-dependent power spectrum from cosmologi¬ 
cal iV-body simulations, and compared the measurements to different theoretical modeling. 
In particular, we have shown that it is not only straightforward to combine the separate 
universe approach with various power spectrum models, but the separate universe approach 
also describes nonlinearity in the squeezed-limit bispectrum due to gravitational evolution 
better than the traditional approach based on the perturbation theory. This would enable 
us to measure the primordial non-Gaussianity in the large-scale structure because we must 
distinguish the primordial signal from the contamination due to the late-time contribution 
that we have computed precisely in this dissertation. 

In chapter [5j we have reported on the first measurement of the three-point function 
with the position-dependent correlation from the SDSS-III BOSS DR10 CMASS sample. 
Since the integrated three-point function of a given subvolume size depends on only one 
separation (unlike the full three-point function which depends on three separations), the 
covariance matrix, which is necessary for the statistical interpretation of the cosmological 
information, can be well estimated with 600 PTHa.los mock catalogs. This allows the 
detection of the amplitude of the three-point function of the BOSS CMASS galaxies at 
7.4(j, which can be turned into the constraint on the nonlinear bias 62 = 0.41 ± 0.41 when 
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combining with the anisotropic clustering and the weak lensing signal. 

We have generalized the study of the response of the small-scale power spectrum to 
m long-wavelength overdensities in chapter [6j This response can be linked to the angle- 
average squeezed-limit (m + 2)-point function, where two modes have wavenumbers much 
larger than the other ones. We have used separate universe simulations, where iV-body 
simulations are performed in the presence of a long-wavelength overdensity by modifying 
the cosmological parameters, to test the separate universe approach on fully nonlinear 
scales to the unprecedented accuracy. We have also tested the standard perturbation 
theory hypothesis that the nonlinear n-point function is completely predicted by the linear 
power spectrum at the equal time. We have found the discrepancies of 10% at k ~ 
0.2 — 0.5 h MpW 1 for five- to three-point functions at z = 0, suggesting that the standard 
perturbation theory fails to describe the correct dynamics of collisionless particles beyond 
these wavenumbers, even if it is calculated to all orders in perturbations. 

While the topic of chapter [6] seems somewhat academic because even measuring the 
three-point function from galaxy surveys is already challenging at the moment, the idea of 
separate universe simulations can largely alleviate the computational resources for studying 
nonlinearities in the squeezed-limit n-point functions. That is, we do not need to perform 
IV-body simulations with a huge volume to simulate the mode coupling between long and 
short wavelength modes. As nonlinearity due to gravitational evolution is the dominant 
contamination for extracting the primordial non-Gaussianity from the large-scale structure 
bispectrum, being able to accurately model the gravity induced bispectrum is currently the 
most important challenge in this field. We can better construct and test the models for 
the squeezed-limit bispectrum with separate universe simulations. 

The quantum origin of all the structures we observe today is one of the most amazing 
ideas in history of physics. Such a bold claim requires careful investigations and validations. 
Upcoming galaxy surveys contain data with unprecedented amount and quality, which 
allow critical tests of this paradigm. The soon-to-be-public BOSS DR12 CMASS sample 
contains approximately 50% more observed galaxies and volume than the DR10 sample. We 
have used the Fisher matrix to show that the BOSS DR12 CMASS sample can potentially 
constrain the local-type primordial non-Gaussianity to be cx(/nl) ~ 17 (95% C.L.). Thus, 
we plan to apply in the near future the same technique to the BOSS DR12 CMASS sample, 
and obtain better constraints on the nonlinear bias, as well as on the logarithmic growth 
rate and the primordial non-Gaussianity. 



Appendix A 

Tree-level redshift-space bispectrum 


In this appendix, we summarize the tree-level redshift-space bispectrum following mm- 


A.l Mapping between real and redshift space 

In redshift space, the observed radial position of an object (galaxy) is the combination of 
the Hubble flow and its peculiar velocity, which is known as redshift-space distortion. The 
mapping between the real-space position x and the redshift-space position s is 

V\\ (x) 

as = ax -|- or s = x — fu\\ (x)£|| , (A.l) 


where x and s are in the comoving coordinates, x\\ is the line-of-sight direction, v\\ — v ■ x\\ 
is the line-of-sight component of the physical peculiar velocity field, u = —v/("H/) is the 
rescaled peculiar velocity field, and T~L = aH = a'/a is the conformal Hubble parameter 
with prime being the derivative with respect to conformal time 


r = 



dt' 

a(t') 


(A.2) 


The density fluctuation in redshift space, <5 s (s), is related to the real-space one, d(x), 
by mass conservation, i.e. 


[1 + 5 s (s)}d 3 s = 

Since the Jacobian, J(x) = d 3 s/d 3 r, is known 
density fluctuation can be written as 


[1 + d(x)]d 3 r . 


(A.3) 


exactly through eq. (A.l), the redshift-space 


_ 1 + J(x) _ J(x) + /V||U||(x) 
J(x) J(x) 


(A.4) 
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where J(x) = 1 — /V||M||(x) and Vy = d/dr\\. In Fourier space, the redshift-space density 
fluctuation is 


d 3 s M±M e - k -[x-/« |( x)i 

J(x) 


s a { k) = I d 3 s 6 S { s)e- iks = 

= [ d 3 r [«5(x) + /V||M||(x)] e ifk Pil( x ) e - ik ' x , 


(A.5) 


where k\\ = k • x\\. Note that in eq. (A.5) the only approximation is the plane-parallel 


approximation, and so it describes the fully nonlinear density fluctuation in redshift space. 
The term in the square brackets describes the so-called “squashing effect”, i.e. the increase 
of the clustering amplitude due to infall into the gravitational potential [80]; the term in 
the exponent encodes the “Finger-of-God effect” which erases power due to the velocity 
dispersion along the line-of-sight [751. 

To proceed, we define the divergence of the rescaled peculiar velocity held as 6*(x) = 
V • u(x), and so 


d 3 r U||(x)e- ik x = 3:11 fl(k) = -^0(k) 


d 3 r V||M||(x)e 


-ik-x 


k 2 

k -h 

k 


k 


W = nle( k), 


(A.6) 


where fj, k = k ■ x\\/k = k\\/k is the cosine of the a ngle b etween k and the line-of-sight. We 
then perturbatively expand e*7 fe ll u lK x ) and use eq. (A.6) to get 


S B (k) = 

= [ d 3 r e“ ,:k x 


d 3 r e * kx [5(x) + /V||U||(x)] \ ^ 
d 3 q 


[ifk ||«||(3 


n =0 


n\ 


x 


i+E 


J (2vr) 3 
iffi k k) n f d 3 q 1 


%) + //#(q) e iq ' x 


n= 1 


(n)! J (27t) 3 


d 3 q n 

(M 


) 0( qi ) • • • ( ) 0(q n )e i(qi+ - +q " ) - x 


[<5(k) +1 Am 

oo 


n=2 


f d 3 q 1 

f d 3 q n r 

"to 

CO 

"to 

CO 


5(qi) + /p^(qi)J [S D \ n 

x wmw—»(<t) • • • —«(q»). 

(n — 1)! q 2 q n 


(A.7) 


where [S D \ n = (27r) 3 h Z) (k - q! 


Qn)- 
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A.2 Redshift-space kernel 

To obtain J s (k), we need to solve the velocity divergence. Since we are interested in scales 
smaller than the Jeans length, we shall treat dark matter and baryons as pressureless fluid. 
Moreover, as the peculiar velocity is much smaller than the speed of light and the scales 
of density fluctuations are much smaller than the horizon size, the system can be treated 
by Newtonian dynamics. The equations of the system are (see, e.g. HZ!) 

S' + V • [(1 + <5)v] = 0 , 
v 7 + (v ■ v) v = -n - V0 , 

v 2 0 = 4 tt Ga 2 pS , (A.8) 


where v = dx/dr = v is the peculiar velocity field in the comoving coordinate with con¬ 
formal time and is equivalent to the physical peculiar velocity field, (j) is the peculiar 
gravitational potential from density fluctuations and p is the mean density. Combining 
eq. (A.8) and the Friedmann equation (i.e. AtiGp(t) = |id 2 fl m (r)), the continuity equa¬ 


tion (the first line in eq. (A.8)) and the Euler equation (the second line in eq. (A.8)) in 
Fourier space can be written as 


d'(k,r) + d(k,r) = 


d 3 h f d 3 ki kk, W1 \ 

[0d] 2—T2-o(k 2 ,T)0(ki,T) 


(2vr) 3 J (2vr) 


d'{ k, r) + U0{ k, r) + -U 2 n m (r)S( k, r) = 


k\ 

d 3 ki 


(2tt) 3 J (2tt) 


d 3 k\ k 2 ( ki-k 2 ) , W1 

[fo]2 0 , 2 , 2 fl(kl,T)fl(k 2 ,T) , 


0 7,2 7,2 

Zj r\j-y r i/2 


(A.9) 


where 6 = V • v is the comoving velocity divergence. Note that 9 = —aHfO. 

To proceed further, we assume that the universe is Einstein de-Sitter, i.e. k2 m (r) 
and a(r) oc r 2 . Then, 5 and 6 can be solved perturbatively as [58> 65] 


= 1 


d(k,r) 

d(k,r) 


OO p 

/ 

n =1 ' 


d 3 qi 
(2vr) 3 ‘ 

a / (r)a n_1 (r) 

77.—1 


f d 3 q n 

J (2vr) 3 

d 3 qi f d 3 q, 


(2tt) ; 


[fe]„Fr ym (qir- - ,qn)^(qi)---^(qn), 

[S D } n Gr ym (qu • • • , qn)^(qi) ■ ■ ■ d z (q„) , 


(2tt) s 


(A.10) 


where Si is the linear density field, iq( insym and G“ nsym are the (unsymmetrized) kernels of 


eq. (A.9), and the recursion relation can be found in [23132!. Strictly speaking, eq. (jAlO]) 
is only valid in Einstein de-Sitter universe, but a good approximation is to replace a = D 
and a' = D 2 Hf (see appendix B.3 in [ 139] ). where D is the linear growth factor and 
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/ = d In D / d In a is the logarithmic growth rate, and eq. (A. 10) can be written as 


9(k,r)=f [ d3qi 1 d3q " 

n= 1 ’’ 


(2vr) 3 J (2 nf 


[&] n C ym ( qi , - ,qn)5i(qi)---5i(qn) , 
[^]nGr ym (qi,---,qnH(qi)---^(qn)- 


(A.ll) 


One can then dehne 

i»(k,r) = D n (r) J |A... J [S D ] n Fr ,m (.qir ■ ■ , 

»«(k,r) =£>"(r) j |A... j J0L [fo]„Gr’“(qi,---,q„)*(qi)---4(q») (A.12) 

such that <5(k, r) = <5 n (k, r) and 9(k, r) = 0„(k, r), and thus <5 n and 9 n are the 

n th order in linear density held. 

For observation, however, one cannot probe dark matter directly, but only the biased 
tracers (e.g. halos or galaxies). Let us assume that the halo (galaxy) density fluctuations 
can be parametrized by the local bias model as |59j 

h(x) = £ V(X) , (A. 13) 

n =0 


where b n are local bias parameters and b 0 assures (Sh) = 0. On the other hand, because of 
the conservation of mass and momentum, we assume that the halo peculiar velocity held 
is identical to the underlying matter peculiar velocity held, i.e. Oh = 9 (see e.g. [26])dJ The 
redshift-space halo density fluctuations can thereby be written as 



d 3 q n 

(2tt) 3 


[fe]n[4(qi) + ft 4#(qi)] 


X 


Udkk) n 1 

(n — 1)! 


q 2 ) ■ ■ ■ 

92 


^9(q„) • 


(A.14) 


■^ote however, that using TV-body simulations |B] recently shows the evidence for linear statistical halo 
velocity bias which remains constant with time, as predicted by the peak model [snug. It is argued 
in mm that the Euler equation has to be modified to correctly describe the coevolution between dark 
matter and halos. 
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Combining eqs. (A.11)-(A.14), one can expand ^, s (k) as 


M k ) =[Mi(k)+/#(k)] 


, x /, x , b 2 f d 3 q 
6i<i 2 (k) + — 


2 J (2vr) 3 
d 3 q 1 f d 3 q 2 


5i(q)<5i(k-q) + /p^ 2 (k) 


[<fo] 2 [Mi(qi) + f t&OiMKf Vkk)— 0i(q 2 


(2tt) 3 ,/ (2tt) 3 


<?2 


+ 


(A.15) 


for which we expand up to the second order of Si. 

For simplicity, one may define the redshift-space kernel Z n ( qi, • • • , q n ) such that 




(2t) 


[$D] n Zn( qi,-" i ( ln)h( C ll)'"b( C l’>) (A. 16 ) 


with 


^i(q 0 

Z 2 (qi,q 2 ) 


bl + f dq, 

biF 2 (qi, q 2 ) + -y + fd 2 G 2 (q 1 , q 2 ) + ^ 




dq 2 


Qi 


<72 


92 


<71' 


(A-17) 


where jl = q • fy/g with q = qi + • • ■ + q n at the n th order, and F n and G n are the 
symmetrized kernels and obtained by taking the mean of all possible permutations of the 
unsymmetrized kernels. Eqs. (A.16)-(A.17) are the main results of this appendix. 

The redshift-space halo bispectrum is defined as 


(^/i,s(ki)^/ liS (k 2 )hh,s(k3)) — (27r) 3 hD(ki + k 2 + k 3 )R/ l)S (ki, k 2 , k 3 ) . (A.18) 

Using the redshift-space kernels, one obtains the leading order terms as 

B h , 8 ( ki, k 2 , k 3 ) = 2[Zi(k 1 )Z 1 (k 2 )Z 2 (ki, k 2 )P l (k 1 )Pi(k 2 ) + 2 cyclic] , (A.19) 


in which we use the fact that Z 2 (— k 3 , —k 2 ) = Z 2 (ki,k 2 ) and Si is Gaussian so that 
(5z(ki)5 z (k 2 )) = (27r) 3 h D (ki + k 2 )Pi(k 1 ) and (S™) = 0 for odd n. 

It is useful to separate B hs into four categories: the linear squashing terms (SQ1) 
which are proportional to P 2 (k l5 k 2 ); the second-order squashing terms (SQ2) which are 
proportional to G 2 (k 1 ,k 2 ); the nonlinear bias terms (NLB) which are related to fe 2 ; the 
damping terms (FOG) due to the velocity dispersion which are not related to P 2 (ki,k 2 ), 
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G 2 { k 1; k 2 ), and b 2 . That is, B hjS = P S qi + P S Q 2 + Pnlb + B FOG such that 


-SsQi(ki, k 2 , k 3 ) 

B SQ 2 (ki, k 2 , k 3 ) 

-Bnlb^i, k 2 , k 3 ) 
-BFOG(ki, k 2 , k 3 ) 


3 

&?^ ( S i - 1 BsQi,i(k 1 ,k2,k 3 ) , 

i =1 
3 

©^^- 1 B SQ2 , i (k 1 ,k 2 ,k 3 ) , 

i =1 
3 

blb 2 ^^/3 l 1 -BNLB,i(ki, k 2 , k 3 ) , 

i —1 

6 3 /5[i?FOG,i(ki, k 2 , k 3 ) + (3B f oG,2(ki, k 2 , k 3 ) + /3PFOG,3(ki, k 2 , k 3 ) 

+ /^-£>FOG,4(ki, k 2 , k 3 ) + /3 2 i?FOG,5(ki, k 2 , k 3 ) + /3 3 -B F oG,6(ki, k 2 , k 3 )] , 

(A.20) 


where /3 = f/b\. The explicit expressions are: for SQ1, 


-B S Qi,i(ki,k 2 ,k 3 ) = 2[F 2 (ki,k 2 )P J (/ci)Pj(/c 2 ) + 2 cyclic] , 

^SQi, 2 (ki, k 2 , k 3 ) = 2[(/i 2 + /i 2 )F 2 (k!, ]t 2 )Pi(ki)Pi(k 2 ) + 2 cyclic] , 

^SQi, 3 (ki, k 2 , k 3 ) = 2[p, 2 p|F 2 (k 1 , k 2 )Pi(k 1 )Pi(k 2 ) + 2 cyclic] ; (A.21) 


for SQ2, 


£sQ2,i(ki,k 2 ,k 3 ) = 2[/r 2 G 2 (k 1 ,k 2 )P i (A;i)P i (/c 2 ) + 2 cyclic] , 

^ S Q2,2(ki, k 2 , k 3 ) = 2[(/i 2 + fil)fi 2 G 2 (ki, k 2 )Pi(ki)Pi(k 2 ) + 2 cyclic] , 

#SQ 2 , 3 (ki, k 2 , k 3 ) = 2[nlnlfi 2 G 2 (k 1 , k 2 )P l (k 1 )Pi(k 2 ) + 2 cyclic] ; (A.22) 


for NLB, 


#NLB,i(ki, k 2 , k 3 ) = Pi(ki)Pi(k 2 ) + 2 cyclic , 
B N LB.2(ki, k 2 , k 3 ) = (/i 2 + nl)Pi(ki)Pi(k 2 ) + 2 cyclic , 
#NLB,3(ki, k 2 , k 3 ) = ^ 2 2 Pi{ki)Pi{k 2 ) + 2 cyclic ; 


(A.23) 
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for FOG, 


-B F oG.i(ki, k 2 , k 3 ) = k/jL ( y + y ) Pi(ki)Pi(k 2 ) + 2 cyclic , 


-BFOG,2(ki, k 2 , k 3 ) —2 


knnin 2 ^ J Pi{ki)Pi(k 2 ) + 2 cyclic 


SFo G) 3 (ki, k 2 , k 3 ) =kfi(^ + P,( W(^) + 2 cyclic , 

^FOG,4(ki, k 2 , k 3 ) = 2 knn\n\ Pi{ki)Pi{k 2 ) + 2 cyclic 

^FOG, 5 (ki,k 2 ,k 3 ) = knnin 2 Pi(ki)Pi(k 2 ) + 2 cyclic , 

-B F oG,6(ki,k 2 ,k 3 ) = knn\nl + jr) p i( k i) p i( k z) + 2 c y clic ■ 


(A.24) 


Note that in eqs. (A.21)-(A.24) we shall simplify the notations // n = Hk n and /i = /j for 
clarity. 


A.3 Tree-level redshift-space integrated bispectrum 
in the squeezed-limit 


Let us define the integrated bispectrum of each component as 

d 2 k f d 3 q a f d 3 q b 


™-mmf 


(2nY 


B x ,i (k - q tt , — k + q a + q b, -q&) 

X W L (q a )W L (-q a - q b )W L (q b ) , (A.25) 


where A" refers to SQ1, SQ2, NLB, or FOG. We numerically evaluate all the components 
of the integrated bispectrum through eq. (A.25) and eqs. (A.21)-(A.24). Figure A.l shows 
the ratios of the components at z = 0 with L = 200 h -1 Mpc. The left and middle 
panels of figure A.l show iBxj(k)/iBx,i(k) for X=SQ1, SQ2, NLB, and FOG with j — 2 
and 3; the right panel shows iBx,i{k)’/iBsQi t i(k) for X=SQ2, NLB, and FOG. We find 
that the ratios become quite scale-independent for k > 0.5 h Mpc -1 , indicating that 


all components have very similar scale-dependencies when the squeezed limit is reached 
(k ^ l/L = 0.005 h Mpc -1 ). Note that in principle the ratios depend on the window 
function and the power spectrum. Fortunately, the anisotropy of the window function can 
be neglected in the squeezed limit. 

To understand the similar scale-dependences for different terms of the tree-level redshift- 
space integrated bispectrum in the squeezed limit, let us now consider the three wavenum¬ 
bers to be ki = k — q a , k 2 = — k + q a + q b, and k 3 = q^ for the ^-configuration of 
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k [h Mpc 1 ] k [h Mpc 1 ] k [h Mpc 1 ] 


Figure A.l: Ratios of the components of the tree-level redshift-space halo bispectrum at 
z — 0 with L = 200 h _1 Mpc, and the line-of-sight direction is z. Left and middle panels: 
iBx,j(k)/iBx,i(k ) for X=SQ1, SQ2, NLB, and FOG, and j — 2 and 3. Right panel: 
iB Xyl (k)/iB SQ1A (k ) for X=SQ2, NLB, and FOG. 


the integrated bispectrum. In the squeezed limit, where k q a ,Qb, we can expand all 
quantities in series of ( q a /k ) and ( qb/k ), and in the leading order (up to ( q a ,b/k)° ) we get 


^1 k ^1 l^ak ^ j ^2 k 1 hafc f-^bk ^ ^ j ^3 % 5 


P l (ki) = P l (k) 


q aJ Uakd 111 Pi(k) 


, Pi(h) = m 


q a [iak + qbi-ibk dlnPi(k) 


k 


din k 


k dink 

Pi(k 3 ) = Pi(q 3 ) , F 2 (k!,k 2 ) = 0 , (j 2 (k 1 ,k 2 ) = 0 , 

77/1 1 \ kqibk q a [^ab . 2 2 I? C\ 1 \ 3 k^L b k <?aha& 2 2 

^(k,.k 3 ) = - 7 -^. fi(k 2 ,k s ) = - +- — -+ . 

/i 1 \ 3 kflbk qa^ab 4 2 /1 1 \ 1 ^h&fc qak'ab 4 2 

G 2 (k 1; k 3 ) = --—-+ ^h bfc , G 2 (k 2 , k 3 ) = — +-—-+ -h fefc , 

Qcihafchfc 9 aha Qahafc T qbk'bk Qahahfc qbk'bk'k 

hi = hfc 4 - 7 - , h 2 = ~Pk H- 7 - , h -3 = -hfe , 


k 


k 


(A.26) 


where hfc = k ■ x \\, ha = <7a • £||, hft = Qb ■ x \\, hfofc = Qb'k and ji ah = q a ■ q b . The redshift-space 
bispectrum can then be expanded in series of ( q a /k ) and (qb/k) as (at the leading order up 
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to ( q a ,b/k)°) 

-®SQl,l,sq 

-BsQl,2,sq 

-SsQl, 3 ,sq 

-BsQ2,l,sq 

-BsQ2,2,sq 

-SsQ 2 , 3 ,sq 

-BNLB,l,sq 

-BNLB,2,sq 

-BNLB, 3 ,sq 
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-BFOG,2,sq 

-BFOG,3,sq 

-BFOG,4,sq 

-®FOG,5,sq 

-®FOG,6,sq 
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jdk Tjdbkdk d~ 2/ifefc/ifc/if, dink ^ bk ^ b ’ 

5 ^ 16 24 ^22 ^ 2 2 2 o 3 

yM/c ^ ^l^bkl^k Y^kl^b “1“ yl Jj bkl J 'kl J/ b ^^bkl^k/^b 
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4 4 din Pl(k) 3 
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a in k 
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where Hk — k ■ xy. Note that eq. (A.27) is independent of q a , so we can simplify the 
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integrated bispectrum as 


iB(k) — >72 / ~n 


d 2 k f d 3 q a f d 3 q b 


V'l J 4tt 


(2?r) 3 J (2vr) 3 


k - q a , —k + q a + q b , -q 6 ) 

X W L (q a )W L (-q a -q b )W L (q b ) 


k^>q a ,qb 


1 f d 2 k f d 3 q, 


= ’ TT 2 / A— / 5 ( k “ Qa, -k + q a + q&, ~qb)W L (q & ) , 


V'l J 47r J (27r) 3 
where we use the fact that 


d 3 g a 

(2tt) 3 


fU L (q a )fU L (-q Q - q 6 ) = W L (q b ) . 


(A.28) 


(A.29) 


If the window function is isotropic, we can angular average both k and q b over eq. (A.27) 
and obtain 


-SsQi,i,sq = Pi(k)Pi(q b ) 

-SsQi,3,sq = Pi(k)Pi(q b ) 
PsQ 2 , 2 ,sq = Pl(k)Pl(q b ) 
PsQ2,3,sq = Pl(k)Pl(q b ) 

-BNLB,2,sq = Pl{k)Pl{q b ) 

-®FOG,l,sq = Pl(k)Pl(q b ) 
-®FOG,3,sq = Pl{k)Pi{q b ) 
-®FOG,5,sq = Pl(k)Pi(q b ) 


47 _ l dlnPf(fc) 

21 3 dlnfc 
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26 dhiPi(k) 
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225 dIn k 
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, -BsQi, 2 , sq — Pi(k)Pi(q b ) 

5 -BsQ 2 ,i,sq = Pi{k)Pi{q b ) 


94 2dlnP;(/c) 
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31 ldlnP/(fc) 
63 9 d In /c 


1225 525 d In k 

A 


, PNLB,l,sq — Pl(k)Pi(q b ) [2] , 


, -B]MLB,3,sq = Pl{k)Pl(q b ) 

J -®FOG,2,sq = Pl(k)Pi(q b ) 
i PFOG,4,sq = Pl(k)Pl(q b ) 

J -®FOG,6,sq = Pl(k)Pl(q b ) 


1 din P;(fc) 

9 d In k 
3 1 dlnP/(fc) 

5 15 d In k 

13 1 d In Pi(k) 

105 _ 21 dlnfc 


22 2 d In P;(fc) 

45 ~~ 15 din A; 

22 2 din P[(k) 

75 ~~ 25 din/c 
13 1 d In P/ (fc) 

175 35 din A; 

(A.30) 


Using eq. (A.30), the integrated bispectrum can further be simplified as 

iB{k) k> ^’ qb PtWalbHxik) , (A.31) 

where af L = yi $ ^I Q b Pi( c lb)Wl(q b ) and Px(k) corresponds to the terms of the four 
categories in the square brackets in eq. (A.30). For n s = 0.95, in the squeezed limit 
~ — 4 = —3.05. Plugging in the value of , we find that the analytical 

(eq. (A.30)) and numerical (figure A.l) results agree well. 























































Appendix B 

Variance of the integrated 
bispectrum estimator 


In this section we compute the variance of the integrated bispectrum estimator. In a 
survey/simulation with volume V r , we first estimate the mean overdensity and the position- 
dependent power spectrum in the j th subvolume with volume Vl by 

h-w- £ (b.i) 

rL 1A reV Lj L kL \lk FL \ek±Ak/2 


where N t l is the number of meshes in V^j, Ar = (Vl/A^l) 1 / 3 is the mesh size, h r .i is the 
discrete density fluctuation held at the integer vector 1, N^l is the number of Fourier modes 
in ( k — A/c/2, k + A/c/2) of Vl , /cfl = 27 t/L is the fundamental frequency of Vl, and Skj 
is the local Fourier transformation of the density fluctuation held in the j th sub-volume. 
The estimated integrated bispectrum is then 

N 3 

-j iy cut 

(B.2) 

JV CUt I 


where N 3 ut = V t /Vl is the number of subvolumes. 

We can rewrite eq. (B.I) using the window function W r3 ,\ {W r] j = 1 if lAr E Vl 3 and 
0 otherwise) as 


5 , = — 
J N, 


rL 


^ ^ S r ,iW r j,i 


lAreVr 


kl 

27r J Vl 


^k,mWkj, m > 


mk F £Vk 


kl 


mA r£VLj 


5 kj ,i = (Ar) 3 Y, ^ ) e in ^ £ 6 k ,j^ n W kj , n , (B.3) 


nk F £V k 


1 /O 

where 14- is the Fourier volume of V r with the fundamental frequency = 27 t/K ' = 
27T/ L r , M r = Ar/cp = 2n/A T r la with = A r / being the mesh number of L r , 
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B. Variance of the integrated bispectrum estimator 


Ml = ArkFL = M r N cut , L = lV cut , and exponential term of 8jk,\ reflects the phase of 
the local Fourier transform (which will then be canceled out when computing the local 
power spectrum Pj(k)). We shall use the capital letters to denote the integer vector mul¬ 
tiplied by N cut . Wki .i is the window function in Fourier space, which can be written as 


Wki,i = (Ar) 3 ^ W r ,, 




= W L (lk F )e- atl " L ‘ , 


(B.4) 


mAr£l4 


where Wl{ k) = Vl IIj= 0 sinc(A;jL/2) is the window function of Vl- Eq. (B.4) is true in 
the continuous limit, i.e. kp —y 0. Namely, there is a slight difference between the dis¬ 
crete Fourier transform of the window function and Wl{ k). However we shall ignore this 
difference and use the continuous limit in the derivation for simplicity. 

Combining the above equations, the integrated bispectrum can be estimated as 


= (| 


VlN tL 


|j|efe±Afe/2 (l,m,n)eVfe 


X 


NL 


E 1 

i =1 


,-i(l+m+n)fc F -r Li 


(B.5) 


Here, to simplify the notation, we drop all the fundamental units (Ar and kp) of the integer 
vectors in the summation and the window function. We find that the p th axis of r Li (r l is 
the center of the z th subvolume) is Tl% v , v — (*p + 1/2)L with i p being the order of the z th 
subvolume in the p th axis, hence the last term of eq. (B.5) can be simplified as 


V 3 ut 

E' 

1=1 


*(l+m+n)fcjr-rij _ 


n<-D 


lp-\-rrip-\-np 


p =0 


sin [(l p + m p + n p ) 7 r] 
sin[(Z p + m p + n p )ir/N cut ] 


(B.6) 


Eq. (B.6) is non-zero only if l p + m p + n p = N cut s p with s p being an integer, and the value 
is (—l) So+Sl+S2 iV c 3 ut for even iV cut and V c 3 ut for odd V cut . We shall assume odd iV cut for 
simplifying the notation, but the results (variance of the estimator) would be identical for 


both cases, as we will show later. We can thus rewrite eq. (B.6) as 


V 3 t) K = /V 3 8 k 

1 'cut^l+m+n.sTVcut 1 v cut u l+m+n,S ) 


(B.7) 


where S^ b is the Kronecker delta and S = sA r cut . Finally, the estimator of the integrated 
bispectrum becomes 


iBL(k)=[ t) VlN kL 


E E E*. 


fc,S-m-n^fc,J-m<5 fe) j + n 


|j|efc±Afc/2 (m,n)eVJt s——oo 


x W L (S - m - n)W L (m)W L (n) . (B.8) 
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Eq. (B.8) is an unbiased estimator because 

9 i 


(iB L (k)) = I ^ 


2tt J V?N, 


1 


L^kL 


E 


E E E 


J+ : 


|j|Efcd=Afc/2 (m,n)GVfc s=—oo 


d 3 q 1 f d 3 q 2 


V?N kL 

L |j|efc±Afc/2' 


(2tt) 3 ,/ (2tt) 3 


x W L (S — m — n)W L ( y m)W L (n) 
B {+qi + q 2 , J k F - q 1; -J k F - q 2 ) 


X W L ( qi )W L (q 2 )W L (- qi -q 2 ) 


N, 


kL 


iB(3k F ) = iB( k) , 


(B.9) 


|j|efc±Afc/2 


where we replace the discrete Fourier transform with the continuous one as [k 3 F /(2'K) 3 ] 
f d 3 q/(2ir) 3 . Note that only s = 0 contributes to ( iB L {k)). 

Replacing the discrete Fourier transform with the integral form, the variance of the 
integrated bispectrum estimator can be computed by 


([W L (k) - (iB L {k))f) = (| \iB L (k )] 1 2 > - [< iB L {k))f 
- (^ F \ 1 V f d 3 q\ C d 3 q 4 

= l 2 n) V r 4 N? r ^ ^ / (27t) 3 J (2vr) 3 

\ / L kL |j 1> j 2 | e A;±Afc/2si,S2=-oo’ / y J J 

(5(ki - qi)5(-ki - q 2 )5(qi + q 2 - gi)<5(k 2 - q 3 )J(-k 2 - q4)5(q3 + q 4 - g 2 )) 

X Rz,(qi)lF L (q2)VF L (gi - qi - q.2)W L (q 3 )W L (q 4 )W L (g 2 - q 3 - Qr) , (B.10) 


where k„ = J n kp and g n = S n k F . We shall assume that the dominant component of the 
six-point function is the disconnected part, and thus the only non-zero component has the 
wavenumber combinations a£] 


k! - qi - k 2 - q 4 = 0 , - ki - q 2 + k 2 - q 3 = 0 , qi + q 2 - gi + q 3 + q 4 - g 2 = 0 . (B.12) 


This gives three delta functions as 


(27r) 9 (5 D (ki - k 2 - q! - q 4 )fe(-ki + k 2 - q 2 - q 2 )5D(gi + g 2 ) • (B.13) 


The last delta function in eq. (B.13) requires gi = —g 2 or s 3 = — s 2 . This means that even 
if even N cut has a parity term (—l) s , it would be canceled because s 4 = —s 2 . 


1 Note that another seemingly non-zero component has the wavenumber combinations as 

ki - qi + k 2 - q 3 = 0 , - ki - q 2 - k 2 - q 4 = 0 , qi + q 2 - gi + q 3 + q 4 - g 2 = 0 . 


(B.ll) 


However, this requires ki = — k 2 or ji = — j 2 . Since we count only the independent Fourier modes, only 

half of the Fourier space is counted (or only the hemisphere is considered), and so it is impossible to have 

jl = ~J2- 
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B. Variance of the integrated bispectrum estimator 


With the above results, the variance of the integrated bispectrum estimator is given by 
(g = S k F = s N cut k F ) 

3 °c 


kf 


2W VtNl L 


E E 

|ji j2|6fc±Afc/2s=-oo ' 


d?q\ f d 3 q 2 


(27r) 3 J (27r) 3 


P( ki - qi)P(-ki - q 2 )P(qi + q 2 - g) 


X W L (qi)W L (q 2 )W L (g - qi - q 2 )Wi(k 1 - k 2 - qi)Vk L (-ki + k 2 - q 2 )VF L (-g + qi + q 2 ) • 

(B.14) 


Note that eq. (B.14) is non-zero only if ^ = k 2 or jx = j 2 , so the variance of the integrated 
bispectrum estimator can be simplified as (k = Jk F = j N cnt k F ) 


1 




E E 

|j|Sfc±Afc/2 s— oo ' 


d 3 q\ f d 3 q 2 


(2t r) 3 J (27t) e 


P(k - qi)P(-k - q 2 )P(qi + q 2 - g) 


x |^( qi )| 2 |W L (q 2 )| 2 |W L (g - qi - q 2 )| 2 , (B.15) 

where we replace k F /(2n) 3 with V r . To proceed further, let us consider the sum over s. 
Since g = s N cut k F = sk FL , we replace the discrete sum with the integral as 

" 3 f d 3 g 


Y p (q-g)|w 7 g-q) 




^F,L 

~2f~ 


^F,L 

27 


(2tt) 


P(q-g)|W L (g-q) 


-3 


T/2 2 _ T/3 

VT <J T — VtO 


L u L 


3 2 
L a L ) 


(B.16) 


where cr 2 = J P (q)\W is the variance of the fluctuation in the volume Vj,. 

Finally, the variance of the integrated bispectrum estimator is simply 


VrV£Nl L 


E *m[/ 


|j|efc±Afc/2 


d 3 q 

(2t r) 3 


V L 


VrNL 


E mwr 


i 3 (k-q)|W 1 (q)| I 


Vi -al\P L (k)f , 


|j|efc±Afc/2 


V r N, 


r^kL 


(B.17) 


where Pl{ k) = dy xhyj P(k—q)|B / z,(q)| 2 is the convolved power spectrum of the subvolume 
Vl- Note that the previous derivation ignores the shot noise contribution. If the shot noise 
is Poisson like, i.e. P s hot = h _1 with n being the number density of the discrete particles, 
then it is trivial to add back as 


Vl 

V r N kL 




[P L (k) + Pshot ] 2 • 


(B.18) 


For the normalized integrated bispectrum, we assume that the variance is dominated by 
the bispectrum term instead of the normalization which contains P^(/c) and cr|, and so its 
variance is 


Vl 

V r N kL 



v L 


cr 


4 

L 


[P L (k ) + Pshot ] 2 

Pl(k) 


(B.19) 
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Figure B.l: Correlation matrices measured from 160 dark matter IV-body simulations 
at z = 0 (left) and 2 (right) for N cut = 4 (bin 0 to 47 in the axis labels), 6 (bin 48 
to 79), 8 (bin 80 to 103), 12 (bin 104 to 119), which correspond to 600, 400, 300, and 
200 h _1 Mpc, respectively. The bottom half and top half of the correlation matrices show 
the unnormalized and normalized integrated bispectrum. 


Figure B.l shows the correlation matrices measured from 160 dark matter IV-body 
simulations at z = 0 and 2. The details of the simulations are given in section |4~T There are 
four subvolume sizes: 600 (bin 0 to 47 in the axis labels), 400 (bin 48 to 79), 300 (bin 80 to 
103), and 200 h Mpc (bin 104 to 119). The bottom (top) half of the correlation matrices 
is the unnormalized (normalized) integrated bispectrum. We find that the cross-correlation 
between different subvolumes is much smaller than that within the same subvolumes. This 
is expected because different subvolumes have different long-wavelength modes, which are 
uncorrelated. At z = 2 the correlation matrices are more diagonal than at z = 0 because 
of the smaller nonlinearity. We also find that the normalization largely diagonalizes the 
correlation matrices, particularly at z = 2 where the off-diagonal components nearly vanish. 
Note that there are stripes, which are more obvious in the normalized integrated bispectrum 
at z — 2, across the diagonal elements of the cross covariances between different sizes of 
subvolumes. This is because these components have the same short-wavelength modes, 
i.e. the scales of the position-dependent power spectrum. Consequently, the correlation is 
stronger. 

Figure [BT2] shows the square root of the variances of the normalized integrated bispectra 
at z = 0 (left) and 2 (right). The data points are measured from simulations, and the solid 
lines are computed using eq. (B.19) with the linear power spectrum and zero shot noise. 
We find that on large scales as well as at high redshift, the agreement between simulations 
and analytical calculation is good. The agreement between the data points and the solid 
lines confirm our calculation of the variance of the integrated bispectrum estimator. 
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B. Variance of the integrated bispectrum estimator 



Figure B.2: Square root of the variances of the normalized integrated bispectrum at z = 0 
(left) and 2 (right). The data points show the measurements from the simulations and the 
solid lines show eq. (B.19) with the linear power spectrum and zero shot noise. 





















Appendix C 

Testing the integrated three-point 
function estimator with Gaussian 
realizations and the local bias model 


We shall demonstrate that our integrated three-point function estimator is unbiased. To 
do this, we first generate the matter density held, S m (r), by Gaussian realization^] with 
P/(k) at z = 0 for the volume V r of 1200 h~ 1 * 3 Mpc 3 and a mesh size of 4 /t _1 Mpc. We then 
compute a mock “halo” density held using the local bias model via 


$h(r) 


Mm(r) + 


h 

2 



E 


rev, ^m(r) 


E 


reVr 


(C.l) 


where we set b\ — 3 and b 2 = 1, and Erey. denotes a sum over grid cells in the entire 
volume. Note that Erev ^( r ) = 0- We then divide the entire volume V r into N s = 12 3 = 
1728 subvolumes Vl of 100 h~ 3 Mpc 3 and N s = 6 3 = 216 subvolumes Vl of 200 h~ 3 Mpc 3 . 
The two-point function in the subvolumes and the integrated three-point function are 
estimated by 


£h(r, r L ) 


Ex+r ,xeV i ^( X + r )^( X ) 


^-E+r,xgVy 


1 


N a 

i(L,h(r) = ^2^ h (r,r L )S h (r L ) , 
i =1 


(C.2) 


where Sh( r L) is the mean halo overdensity in the subvolume centered at r^. Note that 
the denominator in the estimator of £/, (r, r^) takes the boundary effect into account so 
(£h( r , r i)) = €h( r ) without /L,bndry(' r )- This means that the theoretical model of the 
integrated three-point function computed by eq. (2.38) has to be divided by /z^bndryM- 
Since 5 m (r) is Gaussian, the only contribution to the three-point function is from the 
nonlinear bias term, and so the estimated integrated three-point function is exactly given 

<L,b 2 ( r ) 


by 


fL, bndryM ‘ 


1 Since <5 m (r) follows the Gaussian statistics, it is possible that S m ( r) < — 1, which is unphysical. 

However, as we only compute the power spectrum without Poisson sampling the density field, this effect 

can be neglected. 
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C. Testing the integrated three-point function estimator with Gaussian 

realizations and the local bias model 




Figure C.l: The normalized integrated three-point functions of the mock halo density 
field with b\ = 3 and 62 = 1- The left and right panels are for Vl = 100 h ~ 3 Mpc 3 and 
200 h~ 3 Mpc 3 , respectively. The data points show the mean of 300 Gaussian realizations, 
and the error bars are the variances of the mean (but note that the data points are highly 
correlated). The blue dashed and red solid lines are the theoretical models (iCz,,& 2 ) without 
and with [/L,bndry( r )] _1 > respectively. 


Figure C.l shows the normalized integrated three-point functions of the mock halo 
density field with b± = 3 and b -2 — 1 from 300 Gaussian realizations. The measurements are 

<l, 62 _W_ fppjg gives us the confidence that our estimator 


in excellent agreement with 
is unbiased. 


/z.,bndry( r ) ' 











Appendix D 


Effects of effective F2 and G2 kernels 
and non-local tidal bias 


In this appendix, we show how the inferred value of b 2 changes when extending our baseline 
model for the bispectrum based on SPT at the tree level with local bias to the model used 
in the analysis of [62] . 

Their model replaces F 2 and G 2 in eq. (A. 17) with “effective” kernels, (64j and 
G 2 a [63], which are calibrated to match the nonlinear matter bispectrum in of N-body 


simulations. Their model also adds a non-local galaxy bias caused by tidal fields [ 1131 [9, 
EE] to Z 2l i.e., 


Zi — > Z 2 + —b s 


(ki ■ ky 


(D.l) 


where b s 2 = — (4/7)(&i — 1). We use this model to compute the integrated three-point 
function, and find b 2 of the mocks in real and redshift space by performing a joint fit with 


the two-point function as described in section 5.1.3 and 5.1.4 

Table ID. II summarizes the results. The “baseline model” refers to the model based on 
SPT and local bias. The “eff kernel” refers to the model 1 

The “tidal bias” refers to the model with F 2 , G 2l local bias, and tidal bias. Finally, 


Gf, and local bias. 

“both” 


r 2 , local bias, and tidal bias. 


refers to the model with Ff , G| ff , 

Both the effective kernels and the non-local tidal bias result in a larger nonlinear bias, 
which is in better agreement with [62]. The changes of the best-fitting nonlinear bias, 
however, are still within the 1 — a uncertainties, and all the results are consistent with 
[62] . We also calculate the goodness of the fit for all the models in both real and redshift 
space by comparing the mean of the mocks and the best-fitting models, as well as the y 2 - 
distribution. We find that all models perform equally well; thus, we shall primarily use the 
simplest model, i.e. the SPT at the tree level with local bias for modeling the three-point 
function, but also report the results for the extended models. 
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D. Effects of effective F 2 and G 2 kernels and non-local tidal bias 


z-spaee 

bi b 2 

baseline 
eff kernel 
tidal bias 
both 

1.931 ±0.077 0.54 ±0.35 
1.933 ±0.077 0.65 ±0.35 

1.932 ±0.077 0.60 ±0.35 

1.933 ±0.077 0.71 ±0.35 


r-space 

bi b 2 

baseline 
eff kernel 
tidal bias 
both 

1.971 ±0.076 0.58 ±0.31 
1.973 ±0.076 0.62 ±0.31 
1.971 ±0.076 0.64 ±0.31 
1.973 ±0.076 0.68 ±0.31 


Table D.l: Best-fitting values of b\ and b 2 and their uncertainties for mock catalogs, 
obtained using different models of the bispectrum in real space (left) and redshift space 
(right). 





Appendix E 

Comparison for i((r)/a ^ of BOSS 
DR10 CMASS sample and PTHalo 
mock catalogs in different redshift 
bins 


The BOSS DR10 CMASS sample and the mocks have different sets of random samples 
with slightly different n(z ), hence the properties of the observations and the mocks may 
not agree well in all redshift bins. Moreover, as mentioned in [117] . the CMASS sample 
is flux-limited, and thus the observed galaxies statistically have larger stellar masses at 
higher redshift (see figure 1 in [117] ). This may cause redshift evolution of the bias, and so 
the correlation functions. We shall measure iC,{r)/a 2 L as a function of redshift to test this. 


The measurements in the subvolumes are mostly the same as introduced in section 5.1.2 
except that we now measure a(zj) as a function of redshift bin Zj, and the average is done 
in the individual redshift bin. Namely, 


a(zj) = 


E 

E 


iGZj W 9,i 


W, 


r,Zj 


iGz-j 


W , 


9^3 


9fe) = iT-E 




(E.l) 


3 i&z-i 


This assures that S(zj) = 0 for all redshift bins. 

Figure [ eTT] and figure |fh2| show iC,{r)/a 2 L a t different redshift bins for 220 and 120 h _1 Mpc 
subvolumes, respectively. We find no clear sign that iC,{r)/a 2 L of the observations has differ¬ 
ent redshift evolution relative to the mocks. Thus, it is justified to study using one 

effective redshift for the BOSS DR10 CMASS sample. With the upcoming DR12 sample 
with a larger volume, the redshift evolution of fC( r )/u| can be better studied. 
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Figure E.l: *C( r )/ cr i of 220 h 1 Mpc subvolumes in different redshift bins. The redshift 


bins increase from left to right, with the redshift cuts quoted in the beginning of section [5T 
and section 15.21 
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Figure E.2: Same as figure E.l but for 120 h 1 Mpc subvolumes. The redshift bins increase 
from top left to bottom right. 
















































Appendix F 

Squeezed-limit A-point functions and 
power spectrum response 


In this appendix, we prove the relation between the power spectrum response and the 
squeezed limit Appoint functions, as discussed in chapter [6] We only consider equal¬ 
time A-point functions, to which there are no boost-type contributions from kinematic 
consistency relations. Further, we assume that the long-wavelength modes are well inside 
the horizon, removing gauge-dependent terms present for horizon-scale modes. 

We first expand the power spectrum as a function of the linearly extrapolated initial 
overdensity 5 lo as 


P{k,t\S L0 )=Y t ^R n (k,t) \s L0 D(t)T P(k,t) , 
n\ L J 


n =0 


(F.l) 


where R n (k,t) are response functions with Ro(k,t ) = 1. At the same order in derivatives, 
that is at the same order in /q /k of the squeezed-limit A-point function, the power spectrum 
will also depend on the long-wavelength tidal field which can be parametrized through 

KiA k) = f hh - iyi <5(k) . (F.2) 

As here we consider only the angle-averaged long-wavelength modes, the long-wavelength 
tidal field contributions drop out. 

In the following, we will suppress the time argument for clarity. The definition of S n is 
given by 

Sn(k, k'- h, • ■ • , k n ) = J ^ • • • J ^ (5(k)(5(k , )5(k 1 ) • • • S{ k n ))' c , (F.3) 

where the prime denotes that the factor (27r) 3 <5 D (k + k' + k^..^) is dropped and = 
Y^i=i kj. We consider the limit 

S n (k, k'] k ± , ■ ■ ■ , k n ) 

P(k)P l (k 1 )---P l (k n ) ’ 


(F.4) 
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F. Squeezed-limit iV-point functions and power spectrum response 


which means that all |k;| are taken to zero. In this limit, spatial homogeneity enforces 
k' = —k + 0(ki/k), so that (for statistically isotropic initial conditions) the right-hand- 
side of eq. ( F.4[ ) only depends on k. 

In order to prove R n (k ) is equivalent to eq. (F.4), we first note that since we are 
interested in the limit /q —>■ 0, we can replace h(k, ; ) in eq. (F.3) with the linear density held 
5i(kj). We further transform k^ into configuration space, writing 


/ ,l2 h. f 71 [■ 

■■■ j j d 3 Xi e iXi ' kl (h(k)^(k / )h(xi) • • • <5(x n ))' c 

n » 

= \\ J d 3 ' x i e lXi ki S n (k, X\, , x n ) , (F.5) 

i =i ■' 

where 

/ //2 ~ P fj^“T 

v" ■ y • • • 6 ^)c ■ (F-6) 

Note that the angle average is a linear operation and thus commutes with the Fourier 
transform; in other words, the /c-space angle average of the Fourier transform of a function 
is the Fourier transform of the x-space angle average of the same function. 

Now consider the limit hi —>■ 0, which implies that Xi —» oo in the argument of S n . Then 
S n [k ) describes the modulation of the small-scale power spectrum P{k , 0) measured around 
x = 0 by n spherically symmetric large-scale modes (recall that k' « — k). This statement 
can be formalized by introducing an intermediate scale Rl such that 1 /k Rl -C |x, | ~ 
1/ki and defining <5(k) —> h^ L (k) to be the Fourier transform within a cubic volume of size 
Rl around x = 0. Then, ^^(k) = 5(k) + 0(1/(kRi)), while the long-wavelength modes 
are constant over the same volume with corrections suppressed by The corrections 

we expect in the end are thus of order ki/k. To lowest order in these corrections, S n (k ) 
can be written as 


lim : S n (k]x i, • • • ,x n ) = 

ki —^0 


d 2 xi 
47r 


d 2 x n 

47T 


(P(k, 0)5 l (xi) • • • 5l(x„))' c 


(F.7) 


Inserting the expression for the local power spectrum from eq. (F.l), we obtain 


°° 1 f d 2 x i 

lim : S n (k;x i, • • • , x n ) = —R m (k)P(k) ' 


d 2 x 


k-i —^0 


m=0 


m\ 


4n 


4n 


n K(o)s L ( Xl ) ■ ■ ■ s L ( Xn ))'_ 


i- 0 


(R8 > 

Here, the subscript i — 0 indicates that only contractions between 0 and x« are to be taken, 


since the left-hand-side of eq. (F.8) is defined through the connected correlation function 


(all other contractions would contribute to the disconnected part of (h(k)h(k / )h(x 1 ) • • • h(x„))). 


Since all density fields in the correlator in eq. (F.8) are linear, limiting to the contractions 



















F.l Tree-level result: n — 1 
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between 0 and x* then constrains m = n. Therefore, 

1 n 
lirn : S n (k;xi, ■ ■ ■ ,x n ) = —R n (k)P(k)n\Tl 
fci-s-o n\ J -- L 

i =1 


d 2 Xi 

An 




R n (k)P(k)H 


i= 1 


(Pk 

(27T 




(F.9) 


Here, ^ and Pj J denote the linear matter correlation function and power spectrum, respec¬ 
tively. Going back to Fourier space then immediately yields that for k t —> 0, 


S n (k ; Aq, • • • , k n ) = R n {k)P{k) J] P L (ki) + O 


i —1 


h h 


k ’ A; 


NL 


where /cnl is the nonlinear scale. Finally, we can show that 

S n (k, Aq, , kn) 


R n (k ) = lim 


(F.10) 


(F.ll) 


fci-s-o P (A;) Pi,{k\) ■ ■ ■ PL(k n ) 

This provides the connection between the response functions R n {k ) and the angle- 


averaged matter (n + 2)-point function eq. (F.3) in a certain limit squeezed limit (since 
ki -C k). Note that no assumption about the magnitude of k has been made, i.e. this value 
can be fully nonlinear. In the following, we illustrate eq. ( F.ll[ ) at tree level in perturbation 
theory for the cases n = 1 (three-point function) and n = 2 (four-point function). 


F.l Tree-level result: n = 1 


At tree-level for n = 1 we obtain 

lim S^k] ki) tree = evel 2 lim [ 11 [F 2 (k, k x )Pi(k) + F 2 (-k - k 1; k^P^k + ki|)] P;(Aq) . 

fci-»0 fci-s>0 y 47T 

(F.12) 


Eq. (F.ll) then yields 


Ri(k) tree = Ievel —liny [ ^ [P 2 (k, k^P,^) + P 2 (-k - k 1} k^P^k + kx|)] (F.13) 


P/(/c) fci—^0 J 47T 


with 


j-, /, , \ 5 1 (h k 2 \ 2 2 


(F-14) 


where fi is the cosine of kj and k 2 . The term /x/2(A;/Aq) is problematic as we are sending 
ki —> 0. Using that 

h 

k 


|k + ki| = k [1 + qn + 0(q 2 )] , q = — 

din Pi{k) 


Pidk + kil) =P t (k) 


1 + 


din A; 


-'//■I + 0(9 


(F.15) 
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the sum of the two IR-divergent terms in eq. (F.13) becomes 


1 f k id n\ , -( k + k i) ' k i 

2 Rfe PW + Ik + k.lfe 


fl(|k + ki|p = if!(*0 -n 2 


dhiPiik) 

dink 


- 1 


+ O(q) . 
(F.16) 


As expected, the divergent pieces have canceled. We have dropped terms of order k\/k 
which are irrelevant in the limit we are interested in. We finally obtain 


Rl (k) tree = level 2 f 1 dk 
J -i 2 


= 2 


10 

y 


10 1 
y~2 

1 din Pi(k) 
6 d In k 


2 dlnP;(fc) | i 

7 17 I -L 


din k 

1 4 ' 

2 + 21 


+ y 


47 

21 


1 dlnP/(A;) 
3 d In k 


(F-17) 


Eq. (F.17) agrees with the linear prediction for Pi shown in chapter |6j with /i, ei, and g\ 
computed in chapter [3] assuming Einstein-de Sitter universe. 


F.2 Tree-level result: n = 2 

At n = 2, we have 


S n (k;k 1 ,k 2 ) = 

ki,k2—l0 


d 2 k± r d 2 k 2 

47T J 47T 

d 2 ki r d 2 k 2 
47T / 47T 


(5(k)5(k')5(ki)5(k 2 ))' c 

P 3 (k, k 1; k 2 )P z (A:) + P 3 (-k - k 12 , k 1} k 2 )P;(|k + k 12 | 

P 2 (—k 1; k + kx)P 2 (k 2 , k + kx)P;(|k + kx|) 

P 2 (kx, k + k 2 )P 2 (-k 2 , k + k 2 )P(|k + k 2 |) \Pi{ki)Pi(k 2 ) , 

(F.18) 


where P 3 is the symmetrized third-order perturbation theory kernel, and the unsym¬ 
metrized form can be found in [US] . Both F 2 and P 3 contain formally IR-divergent terms 
up to 0[{q^ 2 ) 2 ] where qi >2 = k li2 /k for fcx j2 —* 0 which cancel in the end, so we should 
expand the power spectrum to 0[{ki )2 ) 2 ] to obtain the consistent result at order 0[(gi i2 ) 0 ]. 
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We have 

^(Ik + kil) =P4(A;) 

P(|k + ki + k 2 |) =P(A;) 


9f ri 21 'l k dPi(k) t qln{ k 2 d 2 Pi(k ) | ^^3 


l+(9lA*l + T [l-^l]y P(W dk 


+ 


2 P,(fc) dk 2 

1 + (j9lhl + < 12 ^ 2 } + ^ [9l(l — hi) + 92(1 — hiD 

k dPi(k ) 


+ 0(9? 


29192(^12 - hlh 2 )] 


P/(fc) d In k 


2 (9ihi + 92^2 + 29i92hih2) + 0[(gi, 2 ) 3 ] 


(F.19) 


where /x 1 :2 is the cosine of k and k 12 , and //j 2 is the cosine of k x and k 2 . The leading order 
terms are 


„ /, x tree—level f d k± f C+fc 2 1 
#2(«) = 


47T 


4tt 147 


(628 T 324/ij T 112/1^2 — 280/ii/i 2 /ii 2 

+ 380/1? + 656/i?/i? — 56 /i?/i? 2 ) 
+ (—273/i? + 147/ii/i 2 /ii 2 — 336/i? 


483/i 1 /i 2 + 63/i 2 /t 12 ) 


k dPt(k) 
Pi(k) dk 


147 /i 1 /i 2 


2 2 A; 2 dPP^k) 


P t (k) dk 2 


8420 100 k dPi(k) 1 A; 2 d 2 P L (k ) 

1323 63 p(A;) dk + 9 Pi(k) dk 2 


(F.20) 


Eq. (F.20) also agrees with the linear prediction in chapter [6j with , /*, g t (i = 1,2), 
computed in chapter [3] assuming Eistein-de Sitter universe, are inserted. 
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